%有限差分法的Matlab程序（橢圓型方程） 
function FD_PDE(fun,gun,a,b,c,d)
    % 用有限差分法求解矩形域上的Poisson方程
    tol=10^(-6);  % 誤差界
    N=1000;  % 最大迭代次數
    n=20;  % x軸方向的網格數
    m=20;  % y軸方向的網格數
    h=(b-a)/n; % x軸方向的步長
    l=(d-c)/m; % y軸方向的步長
    for i=1:n-1
        x(i)=a+i*h;
    end % 定義網格點坐標
    for j=1:m-1
        y(j)=c+j*l;
    end % 定義網格點坐標
    u=zeros(n-1,m-1); %對u賦初值
    % 下面定義幾個參數
    r=h^2/l^2;
    s=2*(1+r);
    k=1;
    % 應用Gauss-Seidel法求解差分方程
    while k<=N
        % 對靠近上邊界的網格點進行處理
            % 對左上角的網格點進行處理
            z=(-h^2*fun(x(1),y(m-1))+gun(a,y(m-1))+r*gun(x(1),d)+r*u(1,m-2)+u(2,m-1))/s;
            norm=abs(z-u(1,m-1));
            u(1,m-1)=z;
            % 對靠近上邊界的除第一點和最后點外網格點進行處理
            for i=2:n-2
                z=(-h^2*fun(x(i),y(m-1))+r*gun(x(i),d)+r*u(i,m-2)+u(i+1,m-1)+u(i-1,m-1))/s;
                if abs(u(i,m-1)-z)>norm;
                   norm=abs(u(i,m-1)-z);
                end
                u(i,m-1)=z;
            end
            % 對右上角的網格點進行處理
            z=(-h^2*fun(x(n-1),y(m-1))+gun(b,y(m-1))+r*gun(x(n-1),d)+r*u(n-1,m-2)+u(n-2,m-1))/s;
            if abs(u(n-1,m-1)-z)>norm
               norm=abs(u(n-1,m-1)-z);
            end
            u(n-1,m-1)=z;
        % 對不靠近上下邊界的網格點進行處理
            for j=m-2:-1:2
                % 對靠近左邊界的網格點進行處理
                z=(-h^2*fun(x(1),y(j))+gun(a,y(j))+r*u(1,j+1)+r*u(1,j-1)+u(2,j))/s;
                if abs(u(1,j)-z)>norm
                   norm=abs(u(1,j)-z);
                end
                u(1,j)=z;
                % 對不靠近左右邊界的網格點進行處理
                for i=2:n-2
                    z=(-h^2*fun(x(i),y(j))+u(i-1,j)+r*u(i,j+1)+r*u(i,j-1)+u(i+1,j))/s;
                    if abs(u(i,j)-z)>norm
                       norm=abs(u(i,j)-z);
                    end
                    u(i,j)=z;
                end
                % 對靠近右邊界的網格點進行處理
                z=(-h^2*fun(x(n-1),y(j))+gun(b,y(j))+r*u(n-1,j+1)+r*u(n-1,j-1)+u(n-2,j))/s;
                if abs(u(n-1,j)-z)>norm
                   norm=abs(u(n-1,j)-z);
                end
                u(n-1,j)=z;
            end
        % 對靠近下邊界的網格點進行處理
            % 對左下角的網格點進行處理
            z=(-h^2*fun(x(1),y(1))+gun(a,y(1))+r*gun(x(1),c)+r*u(1,2)+u(2,1))/s;
            if abs(u(1,1)-z)>norm
               norm=abs(u(1,1)-z);
            end
            u(1,1)=z;
            % 對靠近下邊界的除第一點和最后點外網格點進行處理
            for i=2:n-2
               z=(-h^2*fun(x(i),y(1))+r*gun(x(i),c)+r*u(i,2)+u(i+1,1)+u(i-1,1))/s;
               if abs(u(i,1)-z)>norm
                  norm=abs(u(i,1)-z);
               end
               u(i,1)=z;
            end
            % 對右下角的網格點進行處理
            z=(-h^2*fun(x(n-1),y(1))+gun(b,y(1))+r*gun(x(n-1),c)+r*u(n-1,2)+u(n-2,1))/s;
            if abs(u(n-1,1)-z)>norm
               norm=abs(u(n-1,1)-z);
            end
            u(n-1,1)=z;
         % 結果輸出
         if norm<=tol
             fid = fopen('FDresult.txt', 'wt');
             fprintf(fid,'\n********用有限差分法求解矩形域上Poisson方程的輸出結果********\n\n');
             fprintf(fid,'迭代次數: %d次\n\n',k);
             fprintf(fid,'    x的值    y的值       u的值           u的真實值      |u-u(x,y)|\n');
             for i=1:n-1
                 for j=1:m-1
                 fprintf(fid, '%8.3f %8.3f %14.8f   %14.8f  %14.8f\n', [x(i),y(j),u(i,j),gun(x(i),y(j)),abs(u(i,j)-gun(x(i),y(j)))]);
                 end
             end
             fclose(fid);
            break;    % 用來結束while循環
         end
    k=k+1;
    end
    if k==N+1
       fid = fopen('FDresult.txt', 'wt');
       fprintf(fid,'超過最大迭代次數，求解失敗！');
       fclose(fid);