NormalCDF[u_] := CDF[NormalDistribution[0,1],u];
Print[NormalCDF[u]]
Plot[NormalCDF[u],{u,-4,4}]

EuropeanCall[S_,X_,sigma_,r_,d_,tau_]:=Module[{d1,d2},
d1=(Log[S/X]+(r-d+sigma^2/2)*tau)/(sigma*Sqrt[tau]);
d2=d1-sigma*Sqrt[tau];
Return[Exp[-d*tau]*S*NormalCDF[d1]-X*Exp[-r*tau]*NormalCDF[d2]]];

EuropeanPut[S_,X_,sigma_,r_,d_,tau_]:=EuropeanCall[S,X,sigma,r,d,tau]+X*Exp[-r*tau]-S*Exp[-d*tau]

n=5;
Cm=SparseArray[{Band[{1,2}]->-a,Band[{1,1}]->1+2a,Band[{2,1}]->-a},2n-1];
Cm
MatrixForm[Normal[Cm]]

alfa = k/h^2
a = alfa * theta

EuropeanImplicit[u_,h_,k_,n_,m_,theta_]:=Module[{tau,alfa,a,Cm,bm,U},
tau=0.0;
alfa=k/(h^2);
a=alfa*theta;
Cm=SparseArray[{Band[{1,2}]->-a,Band[{1,1}]->1+2a,Band[{2,1}]->-a},2n-1];
For[i=-n,i≤n,i++,x[i]=i*h;u[i,0]=Payoff[x[i],tau]];
For[j=1,j≤m,j++,(*Print["CASOVY KROK CISLO: ",j]*)
tau=j*k;
For[i=-n+1,i<n,i++,
b[i]=u[i,j-1]+alfa*(1.0-theta)*(u[i+1,j-1]-2*u[i,j-1]+u[i-1,j-1])];
u[-n,j]=g[-n]=Lbc[x[-n],tau];
u[n,j]=g[n]=Rbc[x[n],tau];
b[-n+1]+=alfa*theta*g[-n];
b[n-1]+=alfa*theta*g[n];
bm=Table[b[i],{i,-n+1,n-1}];
U=LinearSolve[Cm,bm];
For[i=-n+1,i≤n-1,i++,u[i,j]=U[[i+n]]]
]
]

(*Implicit scheme*)

X=60.0;
r=0.1;
sigma=0.29;
d=0.08;
T=1.0;
k1=r/(0.5*(sigma^2));
k2=(r-d)/(0.5*(sigma^2));
Payoff[x_,tau_]:=Exp[0.25*((((k2-1))^2)+4*k1)*tau]*Max[Exp[0.5*(k2+1)*x]-Exp[0.5*(k2-1)*x],0]
Lbc[x_,tau_]:=0;
Rbc[x_,tau_]:=Exp[0.5*(k2+1)*x+0.25*((((k2+1))^2))*tau]-Exp[0.5*(k2-1)*x+0.25*((((k2-1))^2))*tau]

n=10;
m=5;
h=0.1;
k=T*0.5*(sigma^2)/m;
theta=1.0;
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

g3=ListPlot[Table[{Ssur[i],Vecd[i,0]},{i,-n,n}],PlotStyle->RGBColor[1,0,0],Joined->True];
g4=ListPlot[Table[{Ssur[i],Vecd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g3,g4]

g5(*zelena*)=Plot[EuropeanCall[S,X,sigma,r,d,T],{S,Ssur[-n],Ssur[n]},PlotStyle->RGBColor[0,1,0]];
Show[g3,g4,g5]

g6=ListPlot[Table[{Ssur[i],Vecd[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1]];
Show[g3,g6,g5]

(*ratanie chyby*) Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

ine vstupy

n=20;
m=10;
h=0.05;
k=T*0.5*(sigma^2)/m;
theta=1.0;
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

ine vstupy 2

n=20;
m=20;
h=0.05;
k=T*0.5*(sigma^2)/m;
theta=1.0;
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

Crank nicholson scheme  (theta = 0.5)

n=10;
m=5;
h=0.1;
k=T*0.5*(sigma^2)/m;
theta=0.5; (*TU JE ROZDIEL*)
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

Crank nicholson scheme zjemeny

n=20;
m=10;
h=0.05;
k=T*0.5*(sigma^2)/m;
theta=0.5; (*TU JE ROZDIEL*)
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

0.03812207203328108/0.009582150241250992

(*4 krat lepsie !!! cize radu 2*)

Explicitna schema (theta = 0)

n=10;
m=5;
h=0.1;
k=T*0.5*(sigma^2)/m;
theta=0.0; (*TU JE ROZDIEL*)
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

Explicitna schema zjemnena

n=20;
m=10;
h=0.05;
k=T*0.5*(sigma^2)/m;
theta=0.0; (*TU JE ROZDIEL*)
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

Explicitna schema iny vstup 3

n=10;
m=10;
h=0.1;
k=T*0.5*(sigma^2)/m;
theta=0.0; (*TU JE ROZDIEL*)
EuropeanImplicit[u,h,k,n,m,theta]

g1(*cervena*)=Plot[Payoff[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanCall[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

(*takze teraz je to uuuplne dobre*)

nieco European implicit put

EuropeanImplicitPUT[u_,h_,k_,n_,m_,theta_]:=Module[{tau,alfa,a,Cm,bm,U},
tau=0.0;
alfa=k/(h^2);
a=alfa*theta;
Cm=SparseArray[{Band[{1,2}]->-a,Band[{1,1}]->1+2a,Band[{2,1}]->-a},2n-1];
For[i=-n,i≤n,i++,x[i]=i*h;u[i,0]=PayoffPUT[x[i],tau]];
For[j=1,j≤m,j++,(*Print["CASOVY KROK CISLO: ",j]*)
tau=j*k;
For[i=-n+1,i<n,i++,
b[i]=u[i,j-1]+alfa*(1.0-theta)*(u[i+1,j-1]-2*u[i,j-1]+u[i-1,j-1])];
u[-n,j]=g[-n]=Lbc[x[-n],tau];
u[n,j]=g[n]=Rbc[x[n],tau];
b[-n+1]+=alfa*theta*g[-n];
b[n-1]+=alfa*theta*g[n];
bm=Table[b[i],{i,-n+1,n-1}];
U=LinearSolve[Cm,bm];
For[i=-n+1,i≤n-1,i++,u[i,j]=U[[i+n]]]
]
]

X=60.0;
r=0.1;
sigma=0.29;
d=0.08;
T=1.0;
k1=r/(0.5*(sigma^2));
k2=(r-d)/(0.5*(sigma^2));
PayoffPUT[x_,tau_]:=Exp[0.25*((((k2-1))^2)+4*k1)*tau]*Max[0,-Exp[0.5*(k2+1)*x]+Exp[0.5*(k2-1)*x]]
Rbc[x_,tau_]:=0;
Lbc[x_,tau_]:=-Exp[0.5*(k2+1)*x+0.25*((((k2+1))^2))*tau]+Exp[0.5*(k2-1)*x+0.25*((((k2-1))^2))*tau]

n=10;
m=5;
h=0.1;
k=T*0.5*(sigma^2)/m;
theta=1.0; (*TU JE ROZDIEL*)
EuropeanImplicitPUT[u,h,k,n,m,theta]

g1(*cervena*)=Plot[PayoffPUT[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

(*special thx to dzuraj*)

nieco European put CN

n=10;
m=5;
h=0.1;
k=T*0.5*(sigma^2)/m;
theta=0.5; (*TU JE ROZDIEL*)
EuropeanImplicitPUT[u,h,k,n,m,theta]

g1(*cervena*)=Plot[PayoffPUT[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]

nieco European put CN jemnejsie

n=20;
m=10;
h=0.05;
k=T*0.5*(sigma^2)/m;
theta=0.5; (*TU JE ROZDIEL*)
EuropeanImplicitPUT[u,h,k,n,m,theta]

g1(*cervena*)=Plot[PayoffPUT[x,m*k],{x,-n*h,n*h},PlotStyle->RGBColor[1,0,0]]; 
g2(*modra*)=ListPlot[Table[{x[i],u[i,m]},{i,-n,n}],PlotStyle->RGBColor[0,0,1],Joined->True];
Show[g1,g2]

Table[Table[Vecd[i,j]=X*Exp[-0.5(k2-1)*i*h-(0.25*(((k2-1))^2)+k1)*j*k]*u[i,j],{i,-n,n}],{j,0,m}];
Table[Ssur[i]=X*Exp[x[i]],{i,-n,n}];

Sqrt[Sum[(Vecd[i,m]-EuropeanPut[Ssur[i],X,sigma,r,d,T])^2,{i,-n,n}]/(2n+1)]