源程序：

　　function[u,x,y]=Helmholtz(f,g,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter)

　　%解方程：u_xx+u_yy+g(x,y)u=f(x,y)

　　%自变量取值区域D=[x0,xf,y0,yf]={(x,y)|x0<=x<=xf,y0<=y<=yf}

　　%边界条件

　　%u(x0,y)=bx0(y),u(xf,y)=bxf(y)

　　%u(x,y0)=by0(x),u(x,yf)=byf(x)

　　%x轴均分为Mx段

　　%y轴均分为My段

　　%tol误差因子

　　%MaxIter:最大迭代次数

　　%[u,x,y]:方程u(x，y)在(x,y)点的函数值

　　x0=D(1);xf=D(2);y0=D(3);yf=D(4);

　　dx=(xf-x0)/Mx;x=x0+[0:Mx]*dx;%构造内点数组

　　dy=(yf-y0)/My;y=y0+[0:My]'*dy;

　　Mx1=Mx+1;My1=My+1;

　　%边界条件

　　form=1:My1

　　   u([1Mx1],m)=[bx0(y(m))bxf(y(m))];%左右边界

　　end

　　forn=1:Mx1

　　   u(n,[1My1])=[by0(x(n));byf(x(n))];%上下边界

　　end

　　%边界平均值作迭代初值

　　sum_of_bv=sum(sum([u([1Mx1],2:My)u(2:Mx,[1My1])']));

　　u(2:Mx,2:My)=sum_of_bv/(2*(Mx+My-2));

　　fori=1:Mx

　　   forj=1:My

　　        F(i,j)=f(x(i),y(j));G(i,j)=g(x(i),y(j));

　　   end

　　end

　　dx2=dx*dx;dy2=dy*dy;dxy2=2*(dx2+dy2);

　　rx=dx2/dxy2;ry=dy2/dxy2;rxy=rx*dy2;

　　foritr=1:MaxIter

　　   fori=2:Mx

　　       forj=2:My

　　           u(i,j)=ry*(u(i+1,j)+u(i-1,j))+rx*(u(i,j+1)+u(i,j-1))+rxy*(G(i,j)*u(i,j)-F(i,j));%迭代公式

　　       end

　　   end

　　   ifitr>1&max(max(abs(u-u0)))<MinErr%循环结束条件

　　       break;

　　   end

　　   u0=u;

　　end

　　u=u';

　　在MATLAB中编写脚本文件：

　　f=inline('2*x^2+2*y^2','x','y');

　　g=inline('0','x','y');

　　x0=0;xf=1;y0=0;yf=1;%自变量取值范围

　　Mx=50;My=30;%等分段数

　　bx0=inline('0','y');%边界条件

　　bxf=inline('y^2','y');

　　by0=inline('0','x');

　　byf=inline('x^2','x');

　　D=[x0xfy0yf];MaxIter=100;MinErr=1e-4;

　　[U,x,y]=Helmholtz(f,g,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter);

　　clf,mesh(U)

　　xlabel('x')

　　ylabel('y')

　　zlabel('U')