%function -- blsp_fdcn(X, n, K, r, q, T, m, sigma, lambda)
% European Option pricing using Finite Difference Method
% Crank-Nicolson Scheme
% 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
% m -- the segment number of T
% n -- the segmant number of X
% lambda -- the combination factor

function [x, call, put, Delta, Gamma] ...
    = blsp_fdcn(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;
gamma = 0.5*sigma^2*(x.^2);
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

%{ 
%To compute theta
for i = N:-1:1
    priced = price;
    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
theta = exp(-q*T)*(priced-price)/k;
%}

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