我编写了一个 Matlab 代码，用于在上求解二维泊松方程上具有诺伊曼边界条件， 并且其他边界条件是狄利克雷，我使用了有限差分方法。我的代码如下：$u_{xx} + u_{yy} + f(x,y) = 0$$[a,b]\times [c,d]$$x = b$

function [u,err] = poissonFDM(a,b,c,d,j1,j2,fun,g,typebound,bound,uex,nit)
dltx = (b - a)/(j1+1); dlty = (d - c)/(j2+1);
switch typebound 
    case 'Neumann'
        dim = (j1+1)*j2;
        D = speye(j1+1);
        T = ones(j1+1,1); 
        T = [T,-4*T,T]; T(j1,1) = 2*T(j1,1);
        T = spdiags(T,[-1;0;1],j1+1,j1+1);
        A1 = D;A = T;
        for i = 1 :((j2)- 2)
             A = blkdiag(A,T);
             A1 = blkdiag(A1,D);
        end
        A = blkdiag(A,T);
        A(1:(j1+1)*(j2-1),j1+2:(j1+1)*j2) = A(1:(j1+1)*(j2-1),j1+2:(j1+1)*j2) + A1;
        A(j1+2:(j1+1)*j2,1:(j1+1)*(j2-1)) = A(j1+2:(j1+1)*j2,1:(j1+1)*(j2-1)) + A1;
        h = max(dltx,dlty);
        xh = (a:dltx:b); yh = (c:dlty:d)';
        boundL = g(a,yh);
        boundB = g(xh,c);boundU = g(xh,d);
        gR = bound(b,yh);
        k = 1;
        v = zeros(dim,1);f = zeros(dim,1);Uex = zeros((j1+1)*j2,1);
        for s = 1:j2
             for r = 1:(j1+1)
             f(r,s) = fun(r*dltx + a,s*dlty + c);
                if(nargin == 12)
                    Uex(k) = uex(r*dltx+a,s*dlty+c);
                end
             v(k) = -(h^2)*f(r,s);   
             k = k + 1;
             end
        end
        %Bottom Border     
        v(1:j1) = v(1:j1) - boundB(2:j1+1)';
        v(j1+1) = v(j1+1) - gR(1);
        %Left Border
        v(1:j1:j1*j2 - j1 + 1) = v(1:j1:j1*j2 - j1 + 1) - boundL(2:j2+1);
        %Right Border
        v(j1+1:j1+1:(j1+1)*j2) = v(j1+1:j1+1:(j1+1)*j2) - 2*h*gR(2:j2+1);
        %Up Border
        v((j1+1)*j2-j1:(j1+1)*j2 -1) = v((j1+1)*j2-j1:(j1+1)*j2-1) - boundU(2:j1+1)';
        v((j1+1)*j2)= v((j1+1)*j2 ) - gR(j2+2);
        %solving linear system with conjugate gradiant
        %U = conjgrad(A,v);
        U = A\v;
        %err = U - Uex;
end 
end

我怀疑我的代码的准确性，我想。当我用这个测试我的代码时：$u = \sin(x+y)$

u = poissonFDM(0,pi,0,pi,11,11,@(x,y)2*sin(x+y),@(x,y)sin(x+y),'Neumann',@(x,y)cos(x+y),@(x,y)sin(x+y));
surf(u)

所以图表将是这样的：

的图相去甚远。我的代码有什么问题？关于我使用的方法，我应该说我用中心有限差分公式逼近所有导数！有人可以帮助我吗？$u = sin(x+y)$