In this Matlab code we implement the following finite differences methods: implicit, explicit and Crank-Nicolson to estimate option values. We also plot our concluded values vs. the ones generated by the Black-Scholes closed-form solution.

function FiniteDifferences() S0 = 10; sigma = 0.20;r = 0.04; dT = 0.002;X0 = log(S0);T =0.5 delX = [sigma*sqrt(dT),sigma*sqrt(3*dT),sigma*sqrt(4*dT)]; figure; for i =1:3 dX = delX(i); Path = ceil(log(18/2) / dX) ; Time = (T/dT)+1; X = zeros(Path,1); P = zeros(Path,Time); X(1) = log(2); P(1,Time) = 10 - exp(X(1)); for j =2:Path X(j) = X(j-1) + dX; P(j,Time) = max(10 - exp(X(j)),0); end % Explicit different method [p_d, p_m, p_u] = compute_explicit_prob(r,sigma,dX,dT); A = GetExplicitA(p_u, p_m, p_d,Path); PutValue = zeros(Path,3); % Main Block<span></span> for j = Time-1:-1:1 P(2:Path-1,j) = A*P(:,j+1); P(1,j) = P(2,j) -(exp(X(1)) - exp(X(2))); P(Path,j) = P(Path-1,j) ; end Y = zeros(size(X,1),1); for m =1:size(X,1) [~, Y(m)] = blsprice(exp(X(m)),S0,0.04,0.5,0.20,0); end subplot(3,3,(i-1)*3+1); plot(exp(X),P(:,1),exp(X),Y); legend('Explicit Difference Price','BS Price'); P(:,1:Time-1) = zeros(Path,Time-1); [p_d, p_m, p_u] = compute_implicit_prob(r,sigma,dX,dT); A_I = GetImplicitA(p_u, p_m, p_d,Path); for j = Time-1:-1:1 b = [exp(X(2)) - exp(X(1)); P(2:Path-1,j+1);0 ]; P(:,j) = A_I\b; end subplot(3,3,(i-1)*3+2); plot(exp(X),P(:,1),exp(X),Y); legend('Implicit Difference Price','BS Price'); P(:,1:Time-1) = zeros(Path,Time-1); [p_d, p_m, p_u] = compute_crank_nicol_prob(r,sigma,dX,dT); Z = zeros(Path,1); A_I = GetImplicitA(p_u, p_m, p_d,Path); for j = Time-1:-1:1 Z(2:Path-1) = GetCrankZ(p_u, p_m, p_d,Path,P(:,j+1)); Z(1) = -(exp(X(1)) - exp(X(2))); Z(Path) =Z(Path-1) ; P(:,j) = A_I\Z; end subplot(3,3,(i-1)*3+3); plot(exp(X),P(:,1),exp(X),Y); legend('Crank Nicolson Price','BS Price'); end end function [p_u, p_m, p_d] = compute_explicit_prob(r,sigma,dX,dT) t1 = (sigma^2)/(2*(dX^2)); t2 = (r - ((sigma^2)/2))/(2*dX); t3 = dT * ((sigma^2)/(dX^2)); t4 = r*dT; p_u = dT*(t1 +t2); p_m = 1 - t3 -t4; p_d = dT*(t1 -t2); end function [p_u, p_m, p_d] = compute_implicit_prob(r,sigma,dX,dT) v = r - ((sigma^2)/2); t1 = (sigma^2)/((dX^2)); t2 = (v)/(dX); t3 = dT * ((sigma^2)/(dX^2)); t4 = r*dT; p_u = -0.5* dT*(t1 +t2); p_m = 1 + t3 + t4; p_d = -0.5* dT*(t1 -t2); end function [p_u, p_m, p_d] = compute_crank_nicol_prob(r,sigma,dX,dT) v = r - ((sigma^2)/2); t1 = (sigma^2)/((dX^2)); t2 = (v)/(dX); t3 = dT * ((sigma^2)/(2*dX^2)); t4 = r*dT*0.5; p_u = -0.25* dT*(t1 +t2); p_m = 1 + t3 + t4; p_d = -0.25* dT*(t1 -t2); end function [A] = GetExplicitA(p_u, p_m ,p_d, N) A = zeros(N-2,N); A(1,1:3)= [p_u p_m p_d]; for i = 2:N-2 A(i,:) = circshift(A(i-1,:),[0,1]); end end function [A] = GetImplicitA(p_u, p_m ,p_d, N) A = zeros(N,N); A(1,1:2) = [ 1 -1]; A(2,1:3)= [p_u p_m p_d]; for i = 3:N-1 A(i,:) = circshift(A(i-1,:),[0,1]); end A(N,N-1:N) = [ 1 -1]; end function [A] = GetCrankZ(p_u, p_m ,p_d, N,P) A = zeros(N-2,N); A(1,1:3)= [-p_u -(p_m-2) -p_d]; for i = 2:N-2 A(i,:) = circshift(A(i-1,:),[0,1]); end A = A*P; end