function price = BSCallCN(S0,strike,r,T,vol,h,dt)
% Copyright (c) 2001-2009 by Ales Cerny
% usage BSCallCN(S0,strike,r,T,vol,h,dt)
% computes call option price using Crank-Nicolson finite difference scheme

%************************************************************************%
% BSCallCN.m - supplementary program to                                  %
% Ales Cerny (2009) Mathematical Techniques in Finance (2nd ed.)         % 
% Princeton University Press http://press.princeton.edu/titles/9079.html %
%************************************************************************%
    
% This code is provided 'as-is', without any express or implied warranty.  
% 
% Permission is granted to anyone to use this code for any purpose, 
% subject to the following restrictions:
% 
% 1. The origin of this code must not be misrepresented; you must not
%    claim that you wrote the original code.
% 2. Modified code versions must be plainly marked as such, and must not 
%    be misrepresented as being the original code.
% 3. This notice may not be removed from any source distribution.

% NOTICE TO STUDENTS: To avoid accusations of plagiarism, if you use this 
% code or its modifications in assessed work you should prepend it with a 
% note stating: 
%   "This is the original/modified version of the code BSCallCN.m by 
%    Ales Cerny (2009), Mathematical Techniques in Finance (2nd ed.),  
%    Princeton University Press. The original version is available from 
%    http://www.martingales.info/mtfweb2".
% A similar acknowledgement should appear prominently inside your written 
% report.
  
%***************************%
%       PDE Parameters      %
%***************************%
c = vol*vol;                          % variance of log return          
b = r - c/2;                          % risk-neutral drift of log return

%***************************%
%      grid indexation      %
%***************************%
zeta=8;                               % Number of standard deviations on log stock price grid
Nhalf=ceil(zeta*sqrt(T)*vol/h);       % Number of spatial steps on a positive (negative) half-grid
xgrid=log(S0)+(Nhalf:-1:-Nhalf)*h;
N = 2*Nhalf+1;                        % Number of spatial steps
n = ceil(T/dt);                       % Number of periods
dt=T/n;

%***************************%
%    grid shift to match    %
%       strike value        %
%***************************%
if log(strike)<xgrid(1) && log(strike)>xgrid(N)
    xgrid=xgrid-mod(log(S0/strike),h);
end


%************************%
%    Finite difference   %
%        parameters      %
%************************%     
d1=c*dt/(h*h)/2;
d2=b*dt/h/2;
alph=(d1+d2)/2;
bet=1+d1;
gam=(d1-d2)/2;

if min([alph gam])<0
    disp('Warning BSCallIFD: negative transition probabilities, result will be unreliable');
end

%************************%
%      Precompute LU     %
%      decomposition     %
%************************%    
u = zeros(1,N-2);
u(1) = bet;
for ii = 2 : N-2
    u(ii) = bet - alph*gam/u(ii-1);
end

%************************%
%        payoff          %
%************************%
g = exp(xgrid)-strike;
g = max([g; zeros(1,N)]);

%************************%
%       main loop        %
%************************% 
for tt = 1 : n
    tau=tt*dt;
    gnext=g;                                 
    g(1)=exp(xgrid(1))*exp(r*tau)-strike;
    g(N)=0;
    gtil = alph*gnext(1:(N-2))+(2-bet)*gnext(2:(N-1))+gam*gnext(3:N);
    gtil(1)=gtil(1)+alph*g(1);
    gtil(N-2)=gtil(N-2)+gam*g(N);
    % forward pass solving Ly=gtilde %
    y(1)=gtil(1);
    for ii = 2 : N-3
        y(ii)=gtil(ii)+alph*y(ii-1)/u(ii-1);                    
    end
    ii = N-2;
    y(ii)=gtil(ii)+alph*y(ii-1)/u(ii-1);
    % backward pass solving Ug=y     %
    g(ii+1) = y(ii)/u(ii);
    for ii = N-3 : -1 : 1
        g(ii+1)=(y(ii)+gam*g(ii+2))/u(ii);                    
    end
end 

%************************%
%     interpolation      %
%************************%
price = exp(-r*T)*spline(xgrid,g,log(S0));        % interpolated value using cubic splines %