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

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

NumSteps = tend/h;

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

for i=1:NumSteps,

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

end % for i=1:maxi-1,

figure
plot(T,x)

%CPUTime = etime(clock,bt)


end %The program

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