function V = EuroCall(sigma,r,q,Smax,E,T,M,N)
%
% sigma(i,j) : volatility at (i-1)*delS and t=(j-1)*delt
% delS=Smax/N; delt =T/M;
% r : constant interest rate
% q : continuous dividend yield
% Smax: a large positive underlying price
% T: expiry
% M : number of time steps
% N : number of spatial (price) steps
%
% V: V is the (N+1)-by-(M+1) call value matrix
% V(i,j): put value at S=(i-1)*delS and t=(j-1)*delt
%

V = zeros(N+1,M+1); % initialize price matrix

delS = Smax/N; delt = T/M;

S = linspace(0,Smax,N+1)'; t = linspace(0,T,M+1)'; % grid points

V(1,:) = zeros(1,M+1); % boundary at S =0

V(N+1,:)= Smax*exp(-q*(T-t))'- E*exp(-r*(T-t))'; % boundary at S = Smax

V(:,M+1) = max(S-E,0); % final condition

for j = M:-1:1

e=0.5*(r-q)*delt/delS*S(2:N)-0.5*(delt/delS^2)*sigma(2:N,j).^2.*S(2:N).^2;

d=1+r*delt+delt/delS^2*sigma(2:N,j).^2.*S(2:N).^2;

f=-0.5*(r-q)*delt/delS*S(2:N)-0.5*delt/delS^2*sigma(2:N,j).^2.*S(2:N).^2;

rhs = V(2:N,j+1);
rhs(N-1) = rhs(N-1) - f(N-1)*V(N+1,j);
rhs(1) = rhs(1) - e(1)*V(1,j);

[l,u] = TriDiLU(d,e,f);
temp = LBidiSol(l,rhs);
V(2:N,j) = UBidiSol(u,f,temp);
end;