% matlab code for solving                                                                 %%
% v_tau = v_{yy} + c v_y - d v         (d = 2 r/sigma^2,   c = d-1 for link to options)   %% 
% range:  ymin < y < ymax  .....  take ymin = log(B/K), ymax = ln 3                       %%
%         0 < tau < 1/2 sigma^2 T                                                         %%
% initial cond v(y,0) = (1-e^y)_+                                                         %%
% bdry conds   v(ymin,tau) = 0                                                            %%
%              v(ymax,tau) = 0                                                            %%
% corresponding to a down and out European put option prices with barrier at B            %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r     = 0.05;
sigma = 0.20;
T     = 1;
K     = 10;
d     = 2*r/sigma^2;               % can take any d here
c     = d-1;                       % can take any c here
B     = 5;

ymin  = log(B/K);
ymax  = +log(3);
N     = 100;                       
h     = (ymax-ymin)/N;
k     = h^2/3;                     % initial choice of k - well within stability
M     = ceil(sigma^2*T/2/k);       % "ceil" to make this an integer
k     = sigma^2*T/2/M;             % final choice of k 

y     = ymin:h:ymax;
tau   = 0:k:sigma^2*T/2; 
v     = zeros(N+1,M+1);

v(:,1) = (1 - exp(y)).*(y<0);       % can take any f(y) here 

% solving by the Euler method 
for j=1:M
for i=2:N
    v(i,j+1) = v(i,j) + k/h^2*(v(i+1,j)-2*v(i,j)+v(i-1,j)) + c*k/2/h*(v(i+1,j)-v(i-1,j)) - d*k*v(i,j); 
end
end 

figure(1)    % three-d plot of solution 
Y   = y'*ones(1,M+1);
Tau = ones(N+1,1)*tau;
mesh(Y,Tau,v) 

figure(2)    % three-d plot of associated option prices 
u = K * v;
t = T - 2*Tau/sigma^2; 
x = K * exp(Y);
mesh(x,t,u)

figure(3)    % 2-d plot of option prices for 9 values of tau 
hold on 
plot(x(:,1),u(:,1))
plot(x(:,ceil(M/8)),u(:,ceil(M/8)))
plot(x(:,ceil(2*M/8)),u(:,ceil(2*M/8)))
plot(x(:,ceil(3*M/8)),u(:,ceil(3*M/8)))
plot(x(:,ceil(4*M/8)),u(:,ceil(4*M/8)))
plot(x(:,ceil(5*M/8)),u(:,ceil(5*M/8)))
plot(x(:,ceil(6*M/8)),u(:,ceil(6*M/8)))
plot(x(:,ceil(7*M/8)),u(:,ceil(7*M/8)))
plot(x(:,M),u(:,M))