%Crank-Nicholson格式求解
%求解方程：Ut（x,t）=a1*Uxx(x,t)+a2*Ux(x,t)-r*U(x,t)
%方程的初值及边界条件：U（x，0）=max(exp(x-k,0)
% U(0,t)=0, U(Xmax,t)= V(n)
function [u,x,t] = C_N2(k,xf,T,M,N,z,zz,zzz)
% M ： x 轴的等分段数
% N ： t 轴的等分段数
% s表示方差，r为常数
% x0=3.5 x的最小值
s=0.04;
r=0.04;
dx = xf/M;
x = [0:M]*dx+3.5;
dt = T/N; t = [0:N]'*dt;
for i = 1:M + 1
   u(i,1)=max(exp(x(i))-k,0); %初值条件
end
for n=1:N+1
   D1=(log(exp(x(M+1))/k)+(r+1/2*s)*t(n))/(0.2*t(n)^(1/2));
   D2=D1-(0.2*t(n)^(1/2));
   n1=normcdf(D1,0,1);
   n2=normcdf(D2,0,1);
   V(n)=n1.*(exp(x(M+1)))-k*n2.*(exp(t(n)*(-r)));%边界条件
end
for nn = 1:N +1 %边界设置
    u([1 M + 1],nn) = [0;V(nn)];
end
a1=1/2*s; %系数
a2=r-a1; %系数
r1=(a1/dx/dx-a2/dx)*dt/2; %left
r2 = dt/2*(-2*a1/dx/dx-r) ;%midle
r3=(a1/dx/dx+a2/dx)*dt/2; %right%
R=r3;
e=ones(N+1,1);
P=spdiags([-r1*e (-r2+1)*e -r3*e],[-1 0 1],N+1,M+1); %左矩阵
Q=spdiags([r1*e (r2+1)*e r3*e],[-1 0 1],N+1,M+1); %右矩阵
for K = 2:N + 1
   b=Q(2:M,:)*u(:,K-1)+[zeros(N-2,1);R*V(K)]; % R*V(K)边界条件补充
   u(2:M,K) = P(2:M,2:M)\b;
end
u=u';