我有个偏微分方程的求解问题。这题目应该是用二阶差分方式解得吧。。但是就不知道怎么编程。麻烦您了

function varargout=liu(varargin)
N=5;%改这个N=10 15
a=0;b=1;c=0;d=1;h=1/(N-1);
f=inline('-2*(x^2+y^2)','x','y');
g1x=inline('0');
g2x=inline('x^2');
g1y=inline('0');
g2y=inline('y^2');
chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h);
function [X,Y,Z]=chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h)
%求解下问题
%u_xx+u_yy=-f(x,y) x,y 在区域内 x in a %u=g(x,y) x,y在边界上
%u=g(a,y)=g1y u=g(b,y)=g2y c= %u=g(x,c)=g1x u=g(x,d)=g2x
a %h离散x方向的步长
%h离散y方向的步长
N=10000;
x=a:h:b;
y=c:h:d;
m=length(x);
n=length(y);
ee=0.00001;
[Y,X]=meshgrid(y,x);
Uliu=f0(X,Y);
Z=zeros(m,n);
U=zeros(m,n);
for i=2:m-1
U(i,1)=feval(g1x,x(i));
U(i,n)=feval(g2x,x(i));
end
for j=1:n
U(1,j)=feval(g1y,y(j));
U(m,j)=feval(g2y,y(j));
end
%while true
%下为高斯 赛德尔迭代法
%----------------------------------------------------------------------
for k=1:N

U0=U;

for i=2:m-1

for j=2:n-1

s1=U(i+1,j)+U(i,j+1)+U(i-1,j)+U(i,j-1);

s2=U(i+1,j+1)+U(i-1,j+1)+U(i+1,j-1)+U(i-1,j-1);

s3=feval(f,x(i+1),y(j))+feval(f,x(i-1),y(j))...

+feval(f,x(i),y(j+1))+feval(f,x(i),y(j-1));

U(i,j)=(4*s1+s2-h^2/2*s3-4*h^2*feval(f,x(i),y(j)))/20;

end

end
if max(max(abs(U0-U)))
break
end
end
%-----------------------------------------------------
mesh(X,Y,U)
hold on
plot3(X,Y,Uliu,'ro')
function z=f0(x,y)
%精确解函数
z=(x.*y).^2;