function price = BSCallEFD(S0,strike,r,T,vol,h,dt)
% Copyright (c) 2006-2009 by Ales Cerny
% Usage BSCallEFD(S0,strike,r,T,vol,h,dt)
% Computes call option price using explicit finite difference scheme

%************************************************************************%
% BSCallEFD.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 BSCallEFD.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=5;                               % 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
%dates=0:1:n;
dt=T/n;

%************************%
%    Finite difference   %
%        parameters      %
%************************%     
d1=c*dt/(h*h);
d2=b*dt/h;
pu=(d1+d2)/2;
pm=1-d1;
pd=1-pu-pm; 
if min([pu pm pd])<0
    disp('Warning BSCallEFD: negative transition probabilities, result will be unreliable');
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;
    for ii = 2 : N-1
        g(ii)=pu*gnext(ii-1)+pm*gnext(ii)+pd*gnext(ii+1);                    
    end
end           
price = exp(-r*T)*g(Nhalf+1);