function f=real_fun(x0,t0)

 f=(x0-x0^2)*exp(-t0);

 function f=fun(x0,t0)

 f=(x0^2-x0)*exp(-t0)+2*exp(-t0);

 function f=fi(x0)

 f=x0-x0^2;

 tic;

 clc

 clear

   N=100;

   M=1000;

   t_h=1/M;%t的步长

   x_h=1/N;%x的步长

   x=0:x_h:1;%x的节点

   t=0:t_h:1;%t的节点

  B=-2*ones(1,N-1);

 C=1*ones(1,N-2);

 D=1*ones(1,N-2);

 A=diag(B)+diag(C,1)+diag(D,-1);%三对角矩阵

 F=zeros(N-1,M);

 for i=1:N-1

     for j=1:M

     F(i,j)=fun(x(i+1),t(j));

     end

 end

 F=F.*t_h;

 %********************数值解************************************

 J=-N^2*A*t_h+eye(N-1);%求解线性方程组的系数矩阵

 initial=zeros(N-1,1);

 z=zeros(N-1,M);

 for i=1:N-1

 initial(i)=fi(x(i+1));

 end

 b=initial;

 for j=1:M%控制t的节点

 a=b;

 a=a+F(1:N-1,j);

 z(1:N-1,j)=J\a;%解是n-1维的

 b=z(1:N-1,j);%变成下一次求解的初值

 end

 z=[initial,z];

 Z=[zeros(1,M+1);z;zeros(1,M+1)];

 %********************数值解************************************

 for i=1:N+1

     for j=1:M+1

         real_Z(i,j)=real_fun(x(i),t(j));

     end

 end

 compare=abs(real_Z-Z);

 [X,Y]=meshgrid(x,t);

 % colormap(jet)

 subplot(2,2,1),

 mesh(X,Y,Z');

 xlabel('x');ylabel('t');zlabel('u');title('analytical solution');

 subplot(2,2,2),

 mesh(X,Y,real_Z');

 xlabel('x');ylabel('t');zlabel('u');title('numerical solution');

 subplot(2,2,3),

 mesh(X,Y,compare');

 xlabel('x');ylabel('t');zlabel('u');title('error solution');

 grid on;

 toc;