% Schematic23.m %
clf; clear;start_time=cputime;
sig=0.3; r=0.03; T=1; E=100; dt = 1/360; Nt = T/dt; Smax=300; hstep=2;
a=floor((E-hstep/2)/hstep); b=floor((Smax-E-hstep/2)/hstep);
x=[0 E-(a+0.5)*hstep:hstep:E+(b+0.5)*hstep Smax];
Nx=length(x); h = diff(x); 
h=[h h(end)];u=zeros(Nx,1); f=zeros(1,Nx-1);
main_x(1:Nx-1) = 0; sub_x(1:Nx-1) = 0; sup_x(1:Nx-1) = 0;
lam=2; delta=0.2; k=exp(delta^2/2)-1;
r2=r-lam*k;n=51;xx=linspace(-1,1,n);
G=exp(-0.5*(xx/delta).^2)/(sqrt(2*pi)*delta);
x=[x exp(1)*x(end)+1];
hh=[h(1:end-1) x(end)-x(end-1)];
index(1:Nx,1:n)=-1;
indexmn(1:Nx,1:n)=-1;

for i = 1:Nx
    u(i) = max(x(i)-E,0);
end

for i=2:Nx
    y=x(i)*exp(xx);
    for j=1:n
        index(i,j)=min(find(x>y(j)))-1;
        indexmn(i,j)=(x(index(i,j)+1)-y(j))/hh(index(i,j));
    end
end

for i=2:Nx
    main_x(i-1) = 1/dt+(sig*x(i))^2/(h(i-1)*h(i))...
    -r2*x(i)*(h(i)-h(i-1))/(h(i-1)*h(i))+r;
    sub_x(i-1) = (-(sig*x(i))^2+r2*x(i)*h(i))/(h(i-1)*(h(i-1)+h(i)));
    sup_x(i-1) = -((sig*x(i))^2+r2*x(i)*h(i-1))/(h(i)*(h(i-1)+h(i)));
end

sub_x(Nx-1)=sub_x(Nx-1)-sup_x(Nx-1);
main_x(Nx-1)=main_x(Nx-1)+2*sup_x(Nx-1);

for iter=1:Nt
    u(Nx+1)=u(Nx)+hh(end)/hh(end-1)*(u(Nx)-u(Nx-1));
    for i=2:Nx
        for j=1:n
            z(j)=G(j)*(indexmn(i,j)*u(index(i,j))+(1-indexmn(i,j))*u(index(i,j)+1));
        end
        zG=(xx(2)-xx(1))/3*(z(1)+2*sum(z(3:2:n-2))+4*sum(z(2:2:n-1))+z(n));
        f(i-1)=(1/dt-lam)*u(i)+lam*zG;
    end
    u(2:Nx)=thomas(sub_x,main_x,sup_x,f);
end
end_time=cputime;
time_elapsed=end_time-start_time;
fprintf('Nx = %d \n', Nx)
fprintf('Nt = %d \n', Nt)
fprintf('cputime = %4.8f \n', time_elapsed)
plot(x(1:end-1),u(1:end-1),'k')
axis([0 310 -50 250])
hold on
annotation('line',[0.1297 0.9049],[0.1785 0.1795],'LineWidth',2);
axis off
plot(x(30),u(30),'k')
plot(x(60),u(60),'k')
plot(x(152),u(152),'k')
annotation('line',[0.2722 0.2722],[0.2496 0.1695],'LineStyle',':','LineWidth',1.5);
annotation('line',[0.4228 0.4228],[0.32 0.1695],'LineStyle',':','LineWidth',1.5);
annotation('line',[0.88 0.88],[0.793 0.1695],'LineStyle',':','LineWidth',1.5);
FS=16;
text('Interpreter','latex','String','$${x}/{e}$$','Position',[48 -41 10],'FontSize',FS)
text('Interpreter','latex','String','$$x_i$$','Position',[115 -41 10],'FontSize',FS)
text('Interpreter','latex','String','$$ex$$','Position',[295 -41 10],'FontSize',FS)
plot(x(90),u(90),'ko','MarkerFaceColor','k')
plot(x(100),u(100),'k')
plot(x(110),u(110),'ko')
annotation('line',[0.5727 0.5725],[0.4687 0.1695],'LineStyle',':','LineWidth',1.5);
annotation('line',[0.623 0.623],[0.516 0.1695],'LineStyle',':','LineWidth',1.5);
annotation('line',[0.6723 0.672],[0.567 0.1695],'LineStyle',':','LineWidth',1.5);
text('Interpreter','latex','String','$${x}_k$$','Position',[170 -41 10],'FontSize',FS)
text('Interpreter','latex','String','$$x_ie^{\xi_j}$$','Position',[185 -38 10],'FontSize',FS)
text('Interpreter','latex','String','$$x_{k+1}$$','Position',[210 -41 10],'FontSize',FS)
ss = 10;
i = 30; line([x(i) x(i)],[x(1)-3*ss x(1)-2*ss],'color','k','linewidth',1.5); 
i = 60; line([x(i) x(i)],[x(1)-3*ss x(1)-2*ss],'color','k','linewidth',1.5); 
i = 152; line([x(i) x(i)],[x(1)-3*ss x(1)-2*ss],'color','k','linewidth',1.5); 
i = 90; line([x(i) x(i)],[x(1)-3*ss x(1)-2*ss],'color','k','linewidth',1.5); 
i = 100; line([x(i) x(i)],[x(1)-3*ss x(1)-2*ss],'color','k','linewidth',1.5); 
i = 110; line([x(i) x(i)],[x(1)-3*ss x(1)-2*ss],'color','k','linewidth',1.5); 
i = 90; plot([x(i) x(i)],[x(1)-25 x(1)-25],'ko','MarkerFaceColor','k'); 
i = 110; plot([x(i) x(i)],[x(1)-25 x(1)-25],'ko'); 
text('Interpreter','latex','String','$${u}_k$$','Position',[172 96 10],'FontSize',FS)
text('Interpreter','latex','String','$$u_{ij}$$','Position',[189 111 10],'FontSize',FS)
text('Interpreter','latex','String','$$u_{k+1}$$','Position',[205 132 10],'FontSize',FS)
set(gcf,'PaperPositionMode','auto')