randn('state',100)
tau = 1; 
xi0 = 1;
dt = 0.01;
N = 1000;
T = N*dt;

dW = sqrt(dt)*randn(1,N);         % Brownian increments
W = cumsum(dW);                   % discretized Brownian path 

xi = zeros(1,N);                 % preallocate for efficiency

xi(1) = xi0 - dt*xi0/tau + dW(1); 

for j=2:N
   xi(j) = xi(j-1) - dt*xi(j-1)/tau + dW(j);
end

plot([0:dt:T],[xi0,xi],'b-'), hold off
xlabel('t','FontSize',12)
ylabel('\xi','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
