精确解为：
U（x,t）=e^（x+t）;
此种方法精度为o（h1^2+h2^2）
一：
用追赶法解线性方程组（还可以用迭代法解）
Matlab程序
function[upext]=CN（h1,h2,m,n）
%Crank-Nicolson格式差分法解一维抛物线型偏微分方程
%此程序用的是追赶法解线性方程组
%h1为空间步长，h2为时间步长
%m,n分别为空间，时间网格数
%p为精确解，u为数值解,e为误差
x=（0:
m）*h1+0;x0=（0:
m）*h1;%定义x0,t0是为了f（x,t）~=0的情况%
t=（0:
n）*h2+0;t0=（0:
n）*h2+1/2*h2;
symsf;
for（i=1:
n+1）
for（j=1:
m+1）
f（i,j）=0;%f（i,j）=f（x0（j）,t0（i））==0%
end
end
for（i=1:
n+1）
u（i,1）=exp（t（i））;
u（i,m+1）=exp（1+t（i））;
end
for（i=1:
m+1）
u（1,i）=exp（x（i））;
end
r=h2/（h1*h1）;
for（i=1:
n）%外循环，先固定每一时间层，每一时间层上解一线性方程组%
a
（1）=0;b
（1）=1+r;c
（1）=-r/2;d
（1）=r/2*（u（i+1,1）+u（i,1））+h2*f（i,j）...
+（1-r）*u（i,2）+r/2*u（i,3）;
for（k=2:
m-2）
a（k）=-r/2;b（k）=1+r;c（k）=-r/2;d（k）=h2*f（i,j）+r/2*u（i,k）+（1-r）...
*u（i,k+1）+r/2*u（i,k+2）;
%输入部分系数矩阵，为0的矩阵元素不输入%
end
a（m-1）=-r/2;b（m-1）=1+r;d（m-1）=h2*f（i,j）+r/2*（u（i,m+1）+u（i+1,m+1）...
）+r/2*u（i,m-1）+（1-r）*u（i,m）;
for（k=1:
m-2）%开始解线性方程组消元过程
a（k+1）=-a（k+1）/b（k）;
b（k+1）=b（k+1）+a（k+1）*c（k）;
d（k+1）=d（k+1）+a（k+1）*d（k）;
end
u（i+1,m）=d（m-1）/b（m-1）;%回代过程%
for（k=m-2:
-1:
1）
u（i+1,k+1）=（d（k）-c（k）*u（i+1,k+2））/b（k）;
end
end
for（i=1:
n+1）
for（j=1:
m+1）
p（i,j）=exp（x（j）+t（i））;%p为精确解
e（i,j）=abs（u（i,j）-p（i,j））;%e为误差
end
end
[upext]=CN（0.1,0.005,10,200）;surf（x,t,e）;shadinginterp;
>>xlabel（'x'）;ylabel（'t'）;zlabel（'e'）;
>>title（'误差曲面'）
plot（x,e）
plot（t,e）
误差较向前欧拉法减小一半
但是运行时间较长，约39秒,而前两次运行只需l秒左右；
[upext]=CN（0.01,0.01,100,100）;运行需三分钟左右,误差比前次提高五倍,运算量也提高五倍
[upext]=CN（0.1,0.1,10,10）;surf（x,t,e）运行需要2秒；精度还是挺高的；
[upext]=CN（0.1,0.2,10,5）;surf（x,t,e）
误差还可以接受
此种方法精度高，计算量较大
二：
用迭代法解线性方程组：
Matlab程序如下：
function[uepxtk]=CN1（h1,h2,m,n,kmax,ep）
%解抛物线型一维方程C-N格式（Ut-aUxx=f（x,t）,a>0）
%用g-s（高斯-赛德尔）迭代法解
%kmax为最大迭代次数
%m,n为x,t方向的网格数，例如（2-0）/0.01=200;
%e为误差，p为精确解
symstemp;
u=zeros（n+1,m+1）;
x=0+（0:
m）*h1;
t=0+（0:
n）*h2;
for（i=1:
n+1）
u（i,1）=exp（t（i））;
u（i,m+1）=exp（1+t（i））;
end
for（i=1:
m+1）
u（1,i）=exp（x（i））;
end
for（i=1:
n+1）
for（j=1:
m+1）
f（i,j）=0;
end
end
a=zeros（n,m-1）;
r=h2/（h1*h1）;%此处r=a*h2/（h1*h1）；a=1
for（k=1:
kmax）
for（i=1:
n）
for（j=2:
m）
temp=（（r/2*u（i,j-1）+（1-r）*u（i,j）+r/2*u（i,...
j+1）+h2*f（i,j）+r/2*u（i+1,j-1）+r/2*u（i+1,j+1））/（1+r））;
a（i+1,j）=（temp-u（i+1,j））*（temp-u（i+1,j））;
u（i+1,j）=temp;%此处注意是u（i+1,j）,,而不是u（i+1,j+1）%
end
end
a（i+1,j）=sqrt（a（i+1,j））;
if（k>kmax）
break;
end
if（max（max（a））break;
end
end
for（i=1:
n+1）
for（j=1:
m+1）
p（i,j）=exp（x（j）+t（i））;
e（i,j）=abs（u（i,j）-p（i,j））;
end
end
[uepxtk]=CN1（0.1,0.005,10,200,10000,1e-10）;运行速度：
1秒
迭代次数k=
81
surf（x,t,e）
第二幅图为三角追赶法解方程作出的图，两者几乎一样；
由于迭代法速度很快，所以可以将区间分得更小
[uepxtk]=CN1（0.01,0.01,100,100,10000,1e-12）;surf（x,t,e）;shadinginterp;k=6903
升级会员 