matlab求解微分方程，请大家帮纠错

求解微分方程-Uxx-Uyy+Ux+20yUy+U=f，我已经写出了代码，但迭代第二次误差就很拉垮了，越往后越大，始终检查不出错误。
思路是：fish2d是用来求解-Uxx-Uyy=f类型的快速泊松方程求解器，因此方程化为-Uxx-Uyy=f-Ux-20yUy-U。用初值u0=0代入右端来迭代，每次求解时调用fish2d函数。真解已知，根据真解求出了右端项f。

以下为代码：

1. 主函数

   clc,clear,close all
   a=0;
   b=1;
   c=0;
   d=1;
   N1=31;
   N2=31;
   h1=(b-a)/N1;
   h2=(d-c)/N2;
   x=a+[1:N1-1]*h1;
   y=c+[1:N2-1]*h2;
   [X,Y]=meshgrid(x,y);
   
   %准备矩阵
   e1=ones(N1-1,1);
   e2=ones(N2-1,1);
   e3=ones((N1-1)*(N2-1),1);
   I1=spdiags(e1,0,N1-1,N1-1);
   I2=spdiags(e2,0,N2-1,N2-1);
   I=spdiags(e3,0,(N1-1)*(N2-1),(N1-1)*(N2-1));
   K3=spdiags([-e1,e1],[-1,1],N1-1,N1-1);
   K4=spdiags([-e2,e2],[-1,1],N2-1,N2-1);
   
   %A、B、C、D用以生成系数矩阵
   C=kron(K3,I2);
   C=C/(2*h1);     %u_x系数矩阵
   D=kron(I1,K4);
   D=D/(2*h2);     %u_y系数矩阵
   for i=1:N2-1
       k=i;
       while k<=(N1-1)*(N2-1)
           D(:,k)=D(:,k)*20*y(i);  %a2（x，y）=20y。将矩阵每一列*20*所在列数
           k=k+(N2-1);
       end
   end
   
   %右端项f,代入真解解出
   f=zeros(N1-1,N2-1);
   for i=1:N1-1
       for j=1:N2-1
           f(i,j)=20*y(j)*(10*x(i)*y(j)*exp(x(i)^(9/2))*(x(i)- 1) + 10*x(i)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j) - 1)) - 20*x(i)*exp(x(i)^(9/2))*(x(i)- 1)...
                  - 20*y(j)*exp(x(i)^(9/2))*(y(j) - 1) + 10*x(i)*y(j)*exp(x(i)^(9/2))*(y(j)- 1) + 10*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j)- 1)...
                  - 90*x(i)^(9/2)*y(j)*exp(x(i)^(9/2))*(y(j) - 1) + 10*x(i)*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j) - 1) - (405*x(i)^8*y(j)*exp(x(i)^(9/2))...
                  *(x(i)-1)*(y(j)-1))/2-(495*x(i)^(7/2)*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j)-1))/2+ 45*x(i)^(9/2)*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j) - 1);
       end
   end
   f=reshape(f',(N1-1)*(N2-1),1);%reshape按列取值
   
   % 求解
   u0=zeros((N1-1)*(N2-1),1);
   k=1;
   while k<=2
       u1=fish2d(f-(C+D+I)*u0);
       u0=u1;
       k=k+1;
   end
   %画图
   u1=reshape(u1,N2-1,N1-1);
   mesh(X,Y,u1);
   title('数值解');
   figure
   u_true=10*X.*Y.*(1-X).*(1-Y).*exp(X.^4.5);
   mesh(X,Y,u_true);
   title('真解');
   figure
   mesh(X,Y,u_true-u1);
   title('误差');

2. fish2d函数

   function u = fish2d(f)
   % Poisson solver in 2D based on matlab fft
   % square geometry, homogeneous Dirichlet boundary conditions
   % requires uniform mesh (2^k-1)^2 interior nodes 
   %
   % C. T. Kelley, January 6, 1994
   %
   % This code comes with no guarantee or warranty of any kind.
   %
   % Solves - u_xx - u_yy = f
   %
   % Requires sintv.m and isintv.m
   %
   % function u = fish2d(f)
   %
   %
   global fish_lal fish_ual ua;
   n=length(f);
   m=sqrt(n);
   hh=.5/(m+1);
   h2=(m+1)*(m+1);
   e=ones(m,1);
   %
   % Don't compute eigenvalues, factor tridiagonal matrices
   % or allocate space for the 2-D representation of the solution
   % more than once. Note the use of global variables to store
   % fish_ual, fish_lal, and ua.
   %
   if exist('fish_ual')==1
       [f1, f2]=size(fish_ual);
       m1=sqrt(f1);
   else
       m1=m;
   end
   if (exist('fish_ual') == 0) | (m1 ~= m)
   %
   % allocate room for ua, the 2D representation of the solution
   %
       ua=zeros(m,m);
   %
   % set up diagonal of system using eigenvalues of 1D operator
   %
       x=1:m; x=pi*hh*x'; lam=sin(x); lam=h2*(2*ones(m,1) + 4*lam.*lam);
       en=h2*ones(n,1); fn=en; gn=en; dx=zeros(n,1);
       for i=1:m
           fn(i*m)=0;
           gn((i-1)*m+1)=0;
           dx((i-1)*m+1:i*m)=lam(i)*ones(m,1);
       end
       al=spdiags([-fn, dx, -gn], [-1,0,1], n,n);
   %
   % factor the tridiagonal matrix for solution in the x direction
   %
       [fish_lal,fish_ual]=lu(al);
   end
   %
   % Map f into a rectangular array.
   % Take sine transform in y direction.
   %
   ua(:)=f;
   ua=sintv(ua).';
   %
   % Solve in the x direction
   %
   ua(:)=fish_ual\(fish_lal\ua(:));
   %
   % Take inverse transform.
   %
   ua=isintv((ua.'));
   %
   % Map uz into a linear array.
   %
   u=ua(:);
   
   function u=sintv(z)
   % sine transform, assume that size(z)=[nx,ny]=[2^k-1,ny]
   %
   % C. T. Kelley, May 1994
   %
   % This code comes with no guarantee or warranty of any kind.
   %
   [nx,ny]=size(z);
   ww=-4*fft([zeros(ny,1), z']',2*nx+2);
   u=imag(ww(2:nx+1,:));
   
   function u=isintv(z)
   % inverse sine transform, assume that size(z)=[nx,ny]=[2^k-1,ny]
   % 
   % C. T. Kelley, May 1994
   %
   % This code comes with no guarantee or warranty of any kind.
   %
   [nx,ny]=size(z);
   ww=ifft([zeros(ny,1), z']',2*nx+2);
   u=imag(ww(2:nx+1,:));