function y=BlackScholes_Euler_MC_call(S0,K,T,r,sigma,M,N)
stacksize('max')
    tic();
    Y=rand(M,N,'normal');
    dt=T/N;
    X=zeros(M,N);
    X=cumsum((r-1/2*sigma^2)*dt+sigma*Y*sqrt(dt),2);
    S=S0*exp(X(:,N));
    C=exp(-r*T)*max(S-K,0);
    P=1/M*sum(C);
    VarEst=1/(M-1)*sum((C-P)^2);
    RMSE=sqrt(VarEst)/sqrt(M);
    CI95=[P-1.96*sqrt(VarEst)/sqrt(M),P+1.96*sqrt(VarEst)/sqrt(M)];
    CI99=[P-2.58*sqrt(VarEst)/sqrt(M),P+2.58*sqrt(VarEst)/sqrt(M)];
    time=toc();
    y=[P RMSE; CI95; CI99; time 0]
endfunction