% Black-Scholes simulation
clear
n=10000;
h=1/n;

r=0.1;
sigma=1;
t=linspace(0,30,n+1);

x0=1;
ft(1)=x0;
noice(1)=x0;
Xt(1)=x0;
Bt(1)=0;

dBt=sqrt(h)*randn(1,n);

for i=1:n
    Bt(i+1)=Bt(i)+dBt(i); % Brownian motion
    ft(i+1)=ft(i)*(1+r*(t(i+1)-t(i))); % ODE solution
    noice(i+1)=noice(i)*(1+sigma*dBt(i)); % the noise term
    Xt(i+1)=Xt(i)*(1+r*(t(i+1)-t(i))+sigma*dBt(i)); % the general solution
end

subplot(2,2,1);
plot(t,noice)
subplot(2,2,2);
plot(t,ft)
subplot(2,1,2);
plot(t,Xt)