sims = 100; % number of simulations
dt = 1e-3; % step-size
tv = 0:dt:0.5; % time vector
X0 = 0; % initial condition
No = @(x) normrnd(0,1); % normal random variable
f = @(t,x) sqrt(t)*No(t).*exp(sqrt(t)*No(t)); % function W(t)*exp(W(t))
fmean = @(t) t.*exp(t/2); % Analytic mean

X = zeros(length(tv),1);
X(1) = 0;
Y = zeros(length(X),sims);

% Main Loop
for N = 1:sims
    X = zeros(length(tv),1);
    X(1) = X0;
    for n = 1:length(tv)-1
        X(n+1) = X(n) + (f(tv(n+1),0)-f(tv(n),0));
    end
    Y(:,N) = X;
    disp(N);
end

%% Plots
mean_Y = mean(Y,2);
figure(1);
plot(tv,Y,'g',tv,fmean(tv),'b',tv,mean_Y,'r','LineWidth',2);