function [eul,mil] = q3_24(k,sigma,gam,dt,T)
%compares Euler and Milstein scheme for CEV process, problem 3.24, dt=size of step.
S0=.05;  b=.04;  % initial value =.05 in both cases       
m=T/dt;             % this is the total number of steps
w=cumsum(normrnd(0,sqrt(dt),1,m));   
w1 = [0,w];   % w1 is  a BM on the interval [0,T]
dw=diff(w1);   % these are the increments of the BM
eul = S0*ones(1,m+1);    % initializing
mil = S0*ones(1,m+1);
for i=1:m
    mu=k*(b-eul(1,i)); sig=sigma*(abs(eul(1,i)))^gam;
eul(1,i+1) = eul(1,i) + mu*dt + sig*dw(1,i);
   mu=k*(b-mil(1,i)); sig=sigma*(abs(mil(1,i)))^gam; v=.5*gam*sigma^2*(abs(mil(1,i)))^(2*gam-1);
mil(1,i+1) = mil(1,i)+mu*dt+sig*dw(1,i)+v*((dw(1,i)^2 - dt));
end
x=dt*[0,1:m];
z=.0001*(1:10000);
figure
plot(x,eul,'r',x,mil,'--k')
title(['comparison of Euler (red) and Milstein(dashed) schemes: dt= ' num2str(dt)])