% Euler-Maruyama method for the solution of linear SDEs, with 
% multiplicative noise 

% Wiener part
lambda=2; mu=1; Xzero=1;
T=1; N=2^8; dt=T/N;
dW=sqrt(dt)*randn(1,N);
W=cumsum(dW);

% Exact solution
Xtrue=Xzero*exp((lambda-mu^2/2)*(dt:dt:T)+mu*W);
plot(0:dt:T,[Xzero Xtrue])

% Euler-Maruyama part
R=1; h=R*dt; L=N/R; % L, number of subintervals in EM discretization
Xem=zeros(1,L);
Xtemp=Xzero;
for j=1:L
    Winc=sum(dW((j-1)*R+1:j*R));
    Xtemp=Xtemp+h*lambda*Xtemp+mu*Xtemp*Winc;
    Xem(j)=Xtemp;
end

hold on
plot(0:h:T,[Xzero Xem])