expectedcallmatrix2 = zeros(1,100);
actualcallmatrix2 = zeros(1,100);
imatrix = zeros(1,100);

for i = 1:100  %use a nested loop to obtain different estimated call values
    xvalue2 = zeros(1,10000); %initialize matrix for w values obtained for each iteration at t=.5
    difference2 = zeros(1,10000);  %initialize matrix for difference
    expectedcall2 = zeros(1,10000); %initialize matrix for expected value

    for j = 1:10000 %iterate through at least 10,000 repititions
        [t,w]=milstein([0,.5],12,100);
        xvalue2(j) = w(101); %collect the w value at t=.5 for each iteration
        if xvalue2(j)>15
            difference2(j) = xvalue2(j)-15;
        else
            difference2(j) = 0;
        end
        j = j+1;
    end
    meanxvalue2 = mean(difference2);  
    expectedcallvalue2(i) = exp(-.05*.5)*meanxvalue2;  %compute the expected value of the call
    expectedcallmatrix2(i) = expectedcallvalue2(i);
    actualcallmatrix2(i) = 1;
    imatrix(i) = i;
    i=i+1;
end

actualcallvalue2 = 12*((1+erf(((log(12/15)+(.05+.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)-...  %closed-form solution (fixed) 
   (15*exp(-.05*.5)*(((1+erf(((log(12/15)+(.05-.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)))
actualcallmatrix2 = actualcallvalue2*actualcallmatrix2;

absoluteerror2 = abs(expectedcallmatrix2-actualcallmatrix2);
mm3matrix = mean(expectedcallmatrix2)*ones(1,100);
mm4matrix = mean(absoluteerror2)*ones(1,100);

figure()
plot(imatrix,expectedcallmatrix2,imatrix,actualcallmatrix2,'r',imatrix,mm3matrix,'g')
title('Expected Call Value vs. Iteration Number with the Actual Call Value Superimposed')
xlabel('Iteration Number')
ylabel('Call Value')

figure()
plot(imatrix,absoluteerror2,imatrix,mm4matrix,'r')
title('Absolute Error vs. Iteration Number')
xlabel('Iteration Number')
ylabel('Absolute Error')

errordifference = absoluteerror-absoluteerror2;  
meanerrormatrix = mean(errordifference)*ones(1,100);  

if mean(errordifference) > 0
    'milstein has less error'
else
    'euler-maruyama has less error'
end

figure()
plot(absoluteerror)
hold on
plot(absoluteerror2,'r')
title('Plot of Absolute Errors (Red = Milstein; Blue = Euler-Maruyama')
ylabel('Error')
xlabel('Iteration Number')
hold off

figure()
plot(imatrix,errordifference,imatrix,meanerrormatrix,'r')
title('Relative Error Between Euler-Maruyama and Milstein Methods')
ylabel('Error Difference')
xlabel('Iteration Number')