sigma = 0.4;
T = 1.0;
S0 = 1;
K = 0.8;

N = 2^(0:5);
M = 1000*N^2;

err = zeros(1, length(N));
stat = zeros(1, length(N));
for i=1:length(N) 
    [temp1,temp2] = WeakEulerBS(S0,sigma,K,T,N(i),M(i))
    err(i) = temp1
    stat(i) = temp2
end

coefs = regress(log(1 ./N), log(err))
c = exp(coefs(1))

plot(log(N), log(err), "-", log(N), log(err+stat), ":", log(N), log(err-stat), ":", log(N), log(c./sqrt(N)))