function b=F(f,x0,h,N)

 %       b(1,1)=x0(1)-h*x0(2)-u(1);

 %       b(2,1)=x0(2)+h*x0(1)^2-u(2)-h*f;

 b=zeros(N,1);

 b(1,1)=4*x0(1)-x0(2);

 b(2,1)=h^2*x0(1)^2-2*x0(1)+x0(2)-h^2*f(1)

 for i=3:N

     b(i,1)=x0(i-2)+h^2*x0(i-1)^2-2*x0(i-1)+x0(i)-h^2*f(i-1);

 end

 function g=Jacobian(x0,h,N)

 %   g(1,1)=1;

 %   g(1,2)=-h;

 %   g(2,1)=2*h*x0(1);

 %   g(2,2)=1;

 g=zeros(N);

 g(1,1)=4;

 g(1,2)=-1;

 g(2,1)=2*h^2*x0(1)-2;

 g(2,2)=1;

 for i=3:N

 g(i,i-2)=1;

 g(i,i-1)=2*h^2*x0(i-1)-2;

 g(i,i)=1;

 end

 function x=newton_Iterative_method(f,u,h,N)

  % u：上一节点的数值解或者初值

  % x0 每次迭代的上一节点的数值

  % x1 每次的迭代数值

  % tol 允许误差

  % f 右端函数

 x0=u;

 tol=1e-11;

 x1=x0-Jacobian(x0,h,N)\F(f,x0,h,N);

 while (norm(x1-x0,inf)>tol)

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

     x0=x1;

     x1=x0-Jacobian(x0,h,N)\F(f,x0,h,N);

 end

 x=x1;%不动点

 tic;

 clear

 clc

 N=100;

 h=1/N;

 x=0:h:1;

 y=inline('x.^2.*sin(x).^2+2*cos(x)-x.*sin(x)');

 fun1=inline('x.*sin(x)');

 Accurate=fun1(x);

 f=y(x(2:N));

 Numerical=zeros(N+1,1);

 Numerical(2:end,1)=newton_Iterative_method(f,Numerical(2:end,1),h,N);

 error=Numerical-Accurate';

 subplot(1,3,1)

 plot(x,Accurate);

 xlabel('x');

 ylabel('Accurate');

 grid on

 subplot(1,3,2)

 plot(x,Numerical);

 xlabel('x');

 ylabel('Numerical');

 grid on

 subplot(1,3,3)

 plot(x,error);

 xlabel('x');

 ylabel('error');

 grid on

 toc;