function f=real_fun(x0,t0)

 %精确解

 f=4*x0*(1-x0)*sin(t0);

 function f=F(N,u,U,t,h1,h2)

 %非线性方程组

 %h1是x的步长，h2是t的步长

 %u表示迭代节点,上一时刻的数值解

 %h表示时间节点上的步长

 %N表示空间节点的步数

 a0=0.5*t^4*h2*N^2;

 f(1,1)=a0*(U(2)^2-2*U(1)^2)+h2*fi(h1,t)+u(1)-U(1);

 f(N-1,1)=a0*(-2*U(N-1)^2+U(N-2)^2)+h2*fi((N-1)*h1,t)+u(N-1)-U(N-1);

 for p=2:N-2

 f(p,1)=a0*(U(p+1)^2-2*U(p)^2+U(p-1)^2)+h2*fi(p*h1,t)+u(p)-U(p);

 end

 function f=fi(x0,t0)

 %等式右边的f函数

 f=4*x0*(1-x0)*cos(t0)-16*t0^4*(6*x0^2-6*x0+1)*(sin(t0))^2;

 function g=Jacobian(n,u,t,h1,h2)

 %计算每个时间节点的牛顿迭代过程中的雅可比矩阵

 %u表示迭代初值，上一时刻的数值解作为迭代初值

 a=0.5*t^4*h2*n^2;

 g=zeros(n-1);

 g(1,2)=2*a*u(2);

 g(1,1)=-4*a*u(1);

 g(n-1,n-1)=-4*a*u(n-1);

 g(n-1,n-2)=2*a*u(n-2);

 for p=2:n-2

         g(p,p+1)=2*a*u(p+1);

         g(p,p)=-4*a*u(p);

         g(p,p-1)=2*a*u(p-1);

 end

 g=g-eye(n-1);

 function x=Newtond(n,u,t,h1,h2)

  %使用修改后的牛顿迭代，可以不求雅可比de逆

  %U中间代初值

   %u起始迭代初值

   U=u;

 tol=0.5e-5;

 % Jacobi=Jacobian(n,u,t,h1,h2);%每隔k步求一次雅可比

 x1=U-Jacobian(n,u,t,h1,h2)\F(n,u,U,t,h1,h2);

 while (norm(x1-U,1)>=tol)

     %数值解的1范数是否在误差范围内

     U=x1;

     x1=U-Jacobian(n,u,t,h1,h2)\F(n,u,U,t,h1,h2);

 end

 x=x1;%不动点

 tic;

 clc

 clear

   N=100;

   M=1000;

   t_h=1/M;%t的步长

   x_h=1/N;%x的步长

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

   ti=0:t_h:0.5;%t的节点

 %********************真解**************************

 for i=1:length(x)

     for j=1:length(ti)

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

     end

 end

 %********************真解**************************

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

  ui=zeros(length(x)-2,1);%牛顿迭代初值

 Z=zeros(length(x),length(ti));

 for i=1:length(ti)-1

 Z(2:length(x)-1,i+1)=Newtond(length(x)-1,ui,ti(i+1),x_h,t_h);%t(i+1)时间的牛顿数值解

 ui=Z(2:length(x)-1,i+1);%牛顿迭代初值，上一时刻的数值解作为迭代初值

 end

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

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

 subplot(2,2,1),

 mesh(X,Y,real_Z');

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

 subplot(2,2,2),

 mesh(X,Y,Z');

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

 subplot(2,2,3),

 mesh(X,Y,real_Z'-Z');

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

 title('牛顿迭代法');

 grid on;

 toc;