function y = Black(s0,K,T,sigma)

d1 = log(s0/K)/sigma + sigma/2;

d2 = d1 - sigma;

y = s0 * (0.5 + erf(d1/sqrt(2))/2) - K * (0.5 + erf(d2/sqrt(2))/2);

endfunction