tic;

 clear

 clc

 M=32;%x的步数

 N=16;%y的步数

 h1=1/M;%x的步长

 h2=1/N;%y的步长

 x=0:h1:1;

 y=0:h2:1;

 u=zeros(M+1,N+1);%给数值解分配内存单元

 U=u;%给精确解分配内存单元

 u(1,:)=y.^3;%y边值

 u(M+1,:)=1+y.^3;%y边值

 u(:,1)=x.^3;%x边值。uo

 u(:,N+1)=1+x.^3;%x边值。un

 for i=1:M+1

     for j=1:N+1

     Accurate(i,j)=x(i)^3+y(j)^3;%精确解

     end

 end

 fun=inline('-6*(x+y)','x','y');

 for i=1:M-1

 for j=1:N-1

         f(i,j)=fun(x(i+1),y(j+1));

 end

 end

 Numerical=u;

 error=eye(M+1,N+1);

  while norm(error,inf) >= 1e-5

 for i=2:M

     for j=2:N

         Numerical(i,j)=(f(i-1,j-1)+N^2*u(i,j-1)+M^2*u(i-1,j)+M^2*u(i+1,j)+N^2*u(i,j+1))/(2*M^2+2*N^2);

     end

 end

  error=Numerical-u;

    u=Numerical;

  end

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

 subplot(1,3,1)

 mesh(X,Y,Numerical');

 title('the image of Numerical solution')

 xlabel('x');ylabel('y');zlabel('u');

 subplot(1,3,2)

 mesh(X,Y,Accurate');

 title('the image of Accurate solution')

 xlabel('x');ylabel('y');zlabel('U');

 subplot(1,3,3)

 mesh(X,Y,(Numerical-Accurate)');

 title('the image of error solution')

 xlabel('x');ylabel('y');zlabel('error');

 toc;

 tic;

 clear

 clc

 M=100;%x的步数

 N=100;%y的步数

 h1=1/M;%x的步长

 h2=1/N;%y的步长

 x=0:h1:1;

 y=0:h2:1;

 u=zeros(M+1,N+1);%给数值解分配内存单元

 U=u;%给精确解分配内存单元

 u(1,:)=y.^3;%y边值

 u(M+1,:)=1+y.^3;%y边值

 u(:,1)=x.^3;%x边值。uo

 u(:,N+1)=1+x.^3;%x边值。un

 for i=1:M+1

     for j=1:N+1

     Accurate(i,j)=x(i)^3+y(j)^3;%精确解

     end

 end

 b1=5/3*(M^2+N^2);

 b2=1/6*(5*M^2-N^2);

 b3=1/6*(5*N^2-M^2);

 f=inline('-6*(x+y)','x','y');

 for i=1:M-1

 for j=1:N-1

         ABf(i,j)=(1/144)*(f(x(i),y(j))+10*f(x(i),y(j+1))+f(x(i),y(j+2))...

             +10*f(x(i+1),y(j))+100*f(x(i+1),y(j+1))+10*f(x(i+1),y(j+2))...

             +f(x(i+2),y(j))+10*f(x(i+2),y(j+1))+f(x(i+2),y(j+2)));

 end

 end

 Numerical=u;

 error=eye(M+1,N+1);

  while norm(error,inf) >= 1e-10

 for i=2:M

     for j=2:N

         Numerical(i,j)=(ABf(i-1,j-1)+1/20*b1*Numerical(i-1,j-1)+b3*Numerical(i,j-1)+1/20*b1*Numerical(i+1,j-1)...

             +b2*Numerical(i-1,j)+b2*u(i+1,j)...

             +1/20*b1*u(i-1,j+1)+b3*u(i,j+1)+1/20*b1*u(i+1,j+1))/b1;

     end

 end

  error=Numerical-u;

    u=Numerical;

  end

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

 subplot(1,3,1)

 mesh(X,Y,Numerical');

 title('the image of Numerical solution')

 xlabel('x');ylabel('y');zlabel('u');

 subplot(1,3,2)

 mesh(X,Y,Accurate');

 title('the image of Accurate solution')

 xlabel('x');ylabel('y');zlabel('U');

 subplot(1,3,3)

 mesh(X,Y,(Numerical-Accurate)');

 title('the image of error solution')

 xlabel('x');ylabel('y');zlabel('error');

 toc;