%CH24      Program for Chapter 24
%
%    Crank-Nicolson for a European put 

clf

%%%%%%% Problem and method parameters %%%%%%%
E = 4; sigma = 0.3; r = 0.03; T = 1;
L = 10; Nx = 50; Nt = 50; k = T/Nt; h = L/Nx; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T1 =  diag(ones(Nx-2,1),1) - diag(ones(Nx-2,1),-1);
T2 = -2*eye(Nx-1,Nx-1)  + diag(ones(Nx-2,1),1) + diag(ones(Nx-2,1),-1);
mvec = [1:Nx-1];
D1 = diag(mvec);
D2 = diag(mvec.^2);
F = (1-r*k)*eye(Nx-1,Nx-1) + 0.5*k*sigma^2*D2*T2 + 0.5*k*r*D1*T1;
B = (1+r*k)*eye(Nx-1,Nx-1) - 0.5*k*sigma^2*D2*T2 - 0.5*k*r*D1*T1;
A1 = 0.5*(eye(Nx-1,Nx-1) + F);
A2 = 0.5*(eye(Nx-1,Nx-1) + B);

U = zeros(Nx-1,Nt+1);
U(:,1) = max(E-[h:h:L-h]',0);

for i = 1:Nt
    tau = (i-1)*k;
    p1 = k*(0.5*sigma^2 - 0.5*r)*E*exp(-r*(tau));
    q1 = k*(0.5*sigma^2 - 0.5*r)*E*exp(-r*(tau+k));
    rhs = A1*U(:,i) + [0.5*(p1+q1); zeros(Nx-2,1)];
    X = A2\rhs;
    U(:,i+1) = X;
end

bca = E*exp(-r*[0:k:T]);
bcb = zeros(1,Nt+1);
U = [bca;U;bcb];
mesh([0:k:T],[0:h:L],U)
xlabel('T-t'), ylabel('S'), zlabel('Put Value') 

