%example function call: [t,w]=euler_maruyama([0,.5],12,100);
function[t,w] = milstein(interval,w0,n)
t(1) = interval(1); w(1) = w0;
h = (interval(2)-interval(1))/n;
sqdelt = sqrt(h);
r = 0.05; sigma = 0.25;
for i = 1:n
    db(i) = sqdelt*randn;
    t(i+1) = t(i)+h;
    w(i+1) = w(i)+r*w(i)*h+sigma*w(i)*db(i)+(.5*(sigma^2)*w(i)*((db(i)^2)-h));
end