function [a,b,c,M,N,f,p,u]=finitedifferencemethod2(q,T,S0,K,L,r,sigma,dS, dt) K=100; r=0.02; q=0; sigma=0.2; dS=4; S0=100; T=1; dt=1/200; M=50;%M=L/dS,L=200(nearly double of S0=100); N=200;%N=T/dt,T=1; f=zeros(M+1,N+1); for m=0:M end for m=0:M end for n=0:N end for n=0:N f(1,n+1)=0; end a=zeros(M-1,1); b=a; c=a; for m=1:M-1 f(M+1,n+1)=(M*dS)*exp(-q*(T-n*dt))-K*exp(-r*(T-n*dt)); f(m+1,N+1)=max(m*dS-K,0); f(m+1,1)=S0; end end f(26,1) plot(0:dS:M*dS,f(:,1)) a(m)=1/2*dt*(sigma^2*m^2-(r-q)*m); b(m)=1-dt*(sigma^2*m^2+r); c(m)=1/2*dt*(sigma^2*m^2+(r-q)*m); end for n=N:-1:1 for m=1:M-1 p(m+1,n)=max(0,m*dS-K); u(m+1,n)=a(m)*f(m,n+1)+b(m)*f(m+1,n+1)+c(m)*f(m+2,n+1); f(m+1,n)=max(p(m+1,n),u(m+1,n));