施密特法求解抛物线方程

function varargout=liu(varargin)
%测试问题 来自文献 孙志忠 偏微分数值解法
C=1;a1=0;a2=1;b1=0;b2=1;h1=1/20;h2=1/800;
fx=inline('exp(x)');
gy1=inline('exp(y)');
gy2=inline('exp(1+y)');
[X,Y,U]=Heatflow(fx,gy1,gy2,a1,a2,b1,b2,C,h1,h2);
mesh(X,Y,U);
shading flat;
xlabel('X','FontSize',14);
ylabel('Y','FontSize',14);
zlabel('error','FontSize',14);
title('误差图');
function [X,Y,U]=Heatflow(fx,gy1,gy2,a1,a2,b1,b2,C,h1,h2)
%求解热传导问题：
%CU_y=U_xx
%y(x,b1)=f(x)
%U(a1,y)=g1(y)
%U(a2,y)=g2(y)
%(x,y) in [a1,a2]*[b1,b2]
%m离散x方向的区间个数
%n离散y方向的区间个数
%参考算法：并行算法导论 C.Xavier S.S.lyengar著 张云泉 陈英 译
x=a1:h1:a2;
y=b1:h2:b2;
m=length(x)-1;
n=length(y)-1;
h=(a2-a1)/m;
k=(b2-b1)/n;
r=C*k/h^2;
[X,Y]=meshgrid(x,y);
Z=X;
U=zeros(n+1,m+1);
for i=1:m+1
Z(1,i)=feval(fx,x(i));
end
for i=1:n
Z(i,1)=feval(gy1,y(i));
Z(i,m+1)=feval(gy2,y(i));
end
for i=2:n+1
    for j=2:m
        Z(i,j)=(1-2*r)*Z(i-1,j)+r*(Z(i-1,j-1)+Z(i-1,j+1));
        U(i,j)=abs(Z(i,j)-f(x(j),y(i)));
    end
end
function z=f(x,y)
%精确解函数
z=exp(x+y);