S0 = 1;
K = 1.3;
sigma = 0.1;
r = 0.02;
T = 2; 
N = 24; 
mu = (0:0.05:20)/N ;
a = zeros(2,length(mu)); 
h = T/N; 

M = 50000; 
Z = randn(M,N); 
S = S0 * ones(M,N+1); 

for i=1:length(mu) 
   for j=1:24 
       S(:,j+1) = S(:,j) .* exp( (r-sigma^2/2)*h + sigma*sqrt(h)*(mu(i)+Z(:,j)));
   end
   Smean = mean(S,2); 
   C = exp(-r*T+N*mu(i)^2/2)*(Smean-K).*(Smean>K).*exp(-mu(i)*(N*mu(i)+sum(Z,2))); 
   a(1,i)=100*mean(C);
   a(2,i)=100*std(C)/sqrt(M);
end
   
figure(1)
plot(mu,a(1,:))   
figure(2)
plot(mu,a(2,:))   