function [graph_price,graph_delta,graph_gamma] = ...
    plotfd(X, J, K, r, q, T, N, sigma, lambda, xmin, xmax)

[x, callcn, putcn, Deltacn, Gammacn] ...
    = blsp_fdcn(X, J, K, r, q, T, N, sigma, lambda);
[x, callbdf2, putbdf2, Deltabdf2, Gammabdf2] ...
    = blsp_fdcnbdf2(X, J, K, r, q, T, N, sigma, lambda);
[x, callexp, putexp, Deltaexp, Gammaexp] ...
    = blsp_fdcnexp(X, J, K, r, q, T, N, sigma, lambda);

call = blsprice(x, K, r, T, sigma, q);
delta = blsdelta(x, K, r, T, sigma, q);
gamma = blsgamma(x, K, r, T, sigma, q);

graph_price = figure(1);
plot(x,call,'k', x,callcn,'r', x,callbdf2,'b', x,callexp,'g');
xlim([xmin,xmax]);
xlabel('S, The initial Spot');
ylabel('C(S), The European Call Option Price');
title('Comparation of the 3 different scheme - Price');
legend('The Analytical Solution','Crank-Nicolson Scheme',...
    'Crank-Nicolson Scheme + BDF2',...
    'Crank-Nicolson Scheme + Exponential Fitting');

graph_delta = figure(2);
plot(x,delta,'k', x,Deltacn,'r', x,Deltabdf2,'b', x,Deltaexp,'g');
xlim([xmin,xmax]);
xlabel('S, The initial Spot');
ylabel('Delta(S), The European Call Option Delta');
title('Comparation of the 3 different scheme - Delta');
legend('The Analytical Solution','Crank-Nicolson Scheme',...
    'Crank-Nicolson Scheme + BDF2',...
    'Crank-Nicolson Scheme + Exponential Fitting');

graph_gamma = figure(3);
plot(x,gamma,'k', x,Gammacn,'r', x,Gammabdf2,'b', x,Gammaexp,'g');
xlim([xmin,xmax]);
xlabel('S, The initial Spot');
ylabel('Gamma(S), The European Call Option Gamma');
title('Comparation of the 3 different scheme - Gamma');
legend('The Analytical Solution','Crank-Nicolson Scheme',...
    'Crank-Nicolson Scheme + BDF2',...
    'Crank-Nicolson Scheme + Exponential Fitting');

return