%function -- blsp_fdcnbdf2(X, n, K, r, q, T, m, sigma, lambda)
% European Option pricing using Finite Difference Method
% Crank-Nicolson Scheme + BDF2
% 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_fdcnbdf2(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);

%compute c(N,j)-price1 by 1 step cn method

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;

a1 = matcof*(priceT(1:J-1)');
aJ_1 = L(J-1)*priceX(N+1)+C(J-1)*priceX(N);
a2 = matinv*[zeros(J-2,1);aJ_1];
pricecu = a1 + a2;
price1 = [pricecu; priceX(N)];
price2 = priceT';

%compute c(0,j) by BDF2

B = 3/(2*k) + lambda*2*gammah+r;
H = 2/k - (1-lambda)*2*gammah;

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

for i = N-1:-1:1
    temp = price1;
    p1 = matcof*price1(1:J-1);
    pJ_1 = L(J-1)*priceX(i+1)+C(J-1)*priceX(i);
    pp = [zeros(J-2,1);pJ_1]-0.5/k*price2(1:J-1);
    p2 = matinv*pp;
    pricecut = p1 + p2;
    price1 = [pricecut; priceX(i)];
    price2 = temp;
end

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

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


return