%function -- blsp_fdcnexp(X, J, K, r, q, T, N, sigma, lambda)
% European Option pricing using Finite Difference Method
% Exponential Fitting Scheme + Crank-Nicolson
% Output:
% call -- Call Option price
% put -- Put Option price

% Input:
% K -- Strike
% r -- riskfree rate
% q -- dividen 
% T -- maturity
% sigma -- volatility
% X -- the largest range of x
% N -- the segment number of T
% J -- the segmant number of X
% lambda -- the combination factor

function [x, call, put, Delta, Gamma] ...
    = blsp_fdcnexp(X, J, K, r, q, T, N, sigma, lambda)

h = X/J;
k = T/N;
x = h:h:X;
t = 0:k:T;

ro = r;
r = r-q;
mu = r*x;
gamma1 = 0.5*sigma^2*(x.^2);

if (r == 0 || sigma == 0);
    gamma = gamma1;
elseif (lambda == 0);
    stabtest = tanh(0.5*h*mu./gamma1)./mu;
    if (k/h > min(stabtest));
        error('The method may not be stable. Please use a smaller k/h.');
    end
else
    gamma = 0.5*h*mu.*coth(0.5*h*mu./gamma1);
end


muh = mu/h;
gammah = gamma/(h^2);
A = lambda*(0.5*muh-gammah);
B = 1/k + lambda*2*gammah+r;
C = lambda*(-0.5*muh-gammah);
G = (1-lambda)*(-0.5*muh+gammah);
H = 1/k - (1-lambda)*2*gammah;
L = (1-lambda)*(0.5*muh+gammah);

matrixABC = matrip(A, B, C, J-1);
matrixGHL = matrip(G, H, L, J-1);
matinv = inv(matrixABC);
matcof = matinv*matrixGHL;

priceT = max(0, x-K);
priceX = X - exp(-r*(T-t))*K;
price = priceT';

for i = N:-1:1
    p1 = matcof*price(1:J-1);
    pJ_1 = L(J-1)*priceX(i+1)+C(J-1)*priceX(i);
    p2 = matinv*[zeros(J-2,1);pJ_1];
    pricecut = p1 + p2;
    price = [pricecut; priceX(i)];
end

Delta = diff(price)/h;
Gamma = diff(Delta)/h;
Delta = exp(-q*T)*[Delta;0];
Gamma = exp(-q*T)*[Gamma;0;0];

price = price';
call = exp(-q*T)*price;
forward = exp(-q*T)*x-exp(-ro*T)*K;
put = call - forward;


return