1 tic;
 2 clear
 3 clc
 4 M=[10,20,40,80,160,320,640];%x的步数
 5 K=M; %时间t的步数
 6 for p=1:length(M)
 7 hx=1/M(p);
 8 ht=1/K(p);
 9 r=ht/hx^2; %网格比
10 x=0:hx:1;
11 t=0:ht:1;
12 numerical=zeros(M(p)+1,K(p)+1);
13 numerical(:,1)=exp(x); %初始值
14 numerical(1,:)=exp(t); %边值
15 numerical(M(p)+1,:)=exp(t+1); %边值
16 a=-r/2*ones(M(p)-2,1);b=(1+r)*ones(M(p)-1,1);c=-r/2*ones(M(p)-2,1);
17 fun1=inline('exp(x+t)','x','t');
18 for i=1:length(x)
19     for j=1:length(t)
20 Accurate(i,j)=fun1(x(i),t(j));
21     end
22 end
23 d=r/2*ones(M(p)-2,1);e=(1-r)*ones(M(p)-1,1);f=r/2*ones(M(p)-2,1);
24 B=diag(d,-1)+diag(e,0)+diag(f,1);
25 fun2=inline('0','x','t');
26 for i=1:M(p)-1
27 for k=1:K(p)
28     f(i,k)=ht*fun2(x(i+1),t(k)+ht/2);
29 end
30 end
31 for k=1:K(p)
32 f(1,k)=r/2*(numerical(1,k+1)+numerical(1,k));
33 f(M(p)-1,k)=r/2*(numerical(M(p)+1,k+1)+numerical(M(p)+1,k));
34 end
35 for k=1:K(p)
36     right_vector=f(:,k)+B*numerical(2:M(p),k);
37     numerical(2:M(p),k+1)=chase(a,b,c,right_vector);
38 end
39 error=numerical(2:M(p),2:K(p))'-Accurate(2:M(p),2:K(p))';
40 error_inf(p)=max(max(error));
41 figure(p)
42 [X,Y]=meshgrid(x,t);
43 subplot(1,3,1)
44 mesh(X,Y,Accurate');
45 xlabel('x'),ylabel('t');zlabel('Accurate');
46 title('the image of Accurate result');
47 grid on
48 subplot(1,3,2)
49 mesh(X,Y,numerical');
50 xlabel('x'),ylabel('t');zlabel('numerical');
51 title('the image of numerical result');
52 grid on
53 subplot(1,3,3)
54 mesh(X,Y,numerical'-Accurate');
55 xlabel('x'),ylabel('t');zlabel('error');
56 title('the image of error result');
57 grid on
58 end
59 for k=2:length(M)
60     H=error_inf(p-1)/error_inf(p);
61 E_inf(k-1)=log2(H);
62 end
63 figure(length(M)+1)
64 plot(1:length(M)-1,E_inf,'-r v');
65 ylabel('误差阶数');
66 title('Crank nicolson 误差阶数');
67 grid on
68 toc;