% Schematic14.m %
clear; clf; Smax=300; K=100; r=0.05; sigma=0.2; T=4; delta_tau=0.8; Nt=T/delta_tau;
pw=2; hh=0.5; NS=Smax/hh-mod(Smax,hh)/hh; S=linspace(0,Smax,NS);
for i=1:NS-1
    if (S(i)<=K)&&(S(i+1)>K)
        S(i)=K-hh/2;
        for j=i-1:-1:1
            S(j)=S(j+1)-hh;
        end
        S(1)=0;
        S(i+1)=K+hh/2;
        for j=i+2:NS
            S(j)=S(j-1)+hh;
        end
        S(NS)=Smax;
        break;
    end
end
h=zeros(NS-1,1);
for i=1:NS-1
    h(i)=S(i+1)-S(i);
end
alpha=zeros(NS,1); beta=zeros(NS,1); gamma=zeros(NS,1); f=zeros(NS,1); u=zeros(NS,Nt);
for i=1:NS
    u(i,1)=max(S(i)-K,0)^pw;
end
sum=0; init_S=Smax; N=1000000;
for i=1:N
    S_T=init_S*exp((r-(sigma^2)/2)*T+sigma*randn*sqrt(T));
    sum=sum+exp(-r*T)*(max(S_T-K,0)^pw);
end
Monte=sum/N;
for i=2:NS-1
    alpha(i)=-(((sigma^2)*(S(i)^2))/(h(i-1)*(h(i)+h(i-1))))+((r*S(i)*h(i))/(h(i-1)*(h(i)+h(i-1))));
    beta(i)=(1/delta_tau)+(((sigma^2)*(S(i)^2))/(h(i)*h(i-1)))-((r*S(i)*(h(i)-h(i-1)))/(h(i)*h(i-1)))+r;
    gamma(i)=-(((sigma^2)*(S(i)^2))/(h(i)*(h(i)+h(i-1))))-((r*S(i)*h(i-1))/(h(i)*(h(i)+h(i-1))));
end
for j=1:Nt
    if Monte>0
        u(NS,j+1)=u(NS,1)*exp((1/T)*log(Monte/u(NS,1))*delta_tau*j);
    else
        u(NS,j+1)=0;
    end
    back_up_beta=beta;
    f=u(:,j)/delta_tau;
    f(1)=f(2)-alpha(2)*u(1,j+1);
    f(NS-1)=f(NS-1)-gamma(NS-1)*u(NS,j+1);
    for i=3:NS-1
        beta(i)=beta(i)-((alpha(i)*gamma(i-1))/beta(i-1));
        f(i)=f(i)-((alpha(i)*f(i-1))/beta(i-1));
    end
    u(NS-1,j+1)=f(NS-1)/beta(NS-1);
    for i=NS-2:-1:2
        u(i,j+1)=(f(i)-gamma(i)*u(i+1,j+1))/beta(i);
    end
    beta=back_up_beta;
end
plot(S,u(:,1),'k-','LineWidth',1.3);
axis([0 320 0 80500])
hold on
plot(300,u(NS,:),'ko', 'markerface','k')
text('Interpreter','latex','String','$$\tau=~0~,~u_{N_x}^0$$','Position',[218 40000.0 55.69],'FontSize',18)
text('Interpreter','latex','String','$$\tau=\Delta t,~u_{N_x}^1$$','Position',[218 45800.0 55.69],'FontSize',18)
text('Interpreter','latex','String','$$\vdots$$','Position',[250 63000.0 55.69],'FontSize',20)
text('Interpreter','latex','String','$$\tau=N_t \Delta t,~u_{N_x}^{N_t}=V$$','Position',[183 77000.0 55.69],'FontSize',17)
set(gca, 'fontsize', 18)
set(gcf,'PaperPositionMode','auto')