%
%  SingleEuler.m	
%               Run Euler's method and Milstein on the SDE
%
%                   dX_t = -X_t dt + 2X_t^2 dB_t
%
%  usage:	SingleDiff(tend,initial,h)
%
%               tend = max time for the program to run
%
%               initial = initial condition
%
%               h = size of time discretization for approximate method
%
%
%  example:	SingleDiff(1,1,1/100)
%

function [] = SingleDiff(tend,initial,h)

NumSteps = tend/h;

%initialize.
T = zeros([NumSteps+1,1]);
x_Euler = zeros([NumSteps+1,1]);
x_Milstein = zeros([NumSteps+1,1]);
x_Euler(1) = initial;
x_Milstein(1) = initial;

for i=1:NumSteps,

    R = sqrt(h)*randn;
    x_Euler(i+1) = x_Euler(i) - x_Euler(i)*h + 2*x_Euler(i)^2*R;
    x_Milstein(i+1) = x_Milstein(i) - x_Milstein(i)*h + 2*x_Milstein(i)^2*R + (1/2)*8*x_Milstein(i)^3*(R^2 - h);
    T(i+1) = T(i) + h;

end % for i=1:maxi-1,

figure
plot(T,x_Euler,T,x_Milstein);
legend('Euler-Maruyama','Milstein')

%CPUTime = etime(clock,bt)


end %The program

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
