
n=50
//T=1
delta_t=T/n
//r=0.25
//sigma=0.4
//S0=120
//K=120
up=exp(sigma*sqrt(delta_t))
down=1/up
//R=exp(r*delta_t)
p=(R-down)/(up-down)

function z=cdf(t)
    z=cdfnor("PQ",t,0,1)
endfunction

function y=CallBS(S,K,r,T,sigma)
    d1=(log(S/K)+T*(r+sigma^2/2))/(sigma*sqrt(T))
    d2=d1-sigma*sqrt(T)
    y=S*cdf(d1)-K*exp(-r*T)*cdf(d2)
endfunction


function y=PutBS(S,K,r,T,sigma)
    d1=(log(S/K)+T*(r+sigma^2/2))/(sigma*sqrt(T))
    d2=d1-sigma*sqrt(T)
    y=K*exp(-r*T)*cdf(-d2)-S*cdf(-d1)
endfunction
