如有问题联系:1759553924
function [err,zuidaerr,U,Ujing]=houeulererwei(m,n,jingqueu,phi,f)
%此程序用于求解差分方程du/dt-d2u(x)-d2u(y)=f(x)的解
%u(x,0)=phi,u(a,t)=alpha,u(b,t)=beta
%差分空间剖分m等分,时间剖分n等分

h=(1-0)/m;
t=(1-0)/n;
row=t/h;
r=t/h.^2;
x=0:h:1;
alpha=0.001*ones(1,m+1);
y=0:t:1;
for i=1:m+1
U(i,1:m+1)=feval(phi,(i-1)*h,x);
end
U(1:m+1,1)=zeros(m+1,1);
U(1:m+1,1+m)=zeros(m+1,1);
U(1,1:m+1)=zeros(1,1+m);
U(m+1,1:1+m)=zeros(1,1+m);

   A=(1/h/2)*ones(1,m-2)-[alpha(2:m-1)+alpha(3:m)]/2/h^2;
   C=-1/h/2*ones(1,m-2)-[alpha(2:m-1)+alpha(3:m)]/2/h^2;
   B=(1/t)*ones(1,m-1)+[2*alpha(2:m)+alpha(1:m-1)+alpha(3:m+1)]/h^2;
 AA=diag(B)+diag(C,-1)+diag(A,1);
 AAAA=[];
 for k=1:m-1
     AAAA=blkdiag(AAAA,AA);
 end
   A1=(1/h/2)*ones(1,m-1)-[alpha(2:m)+alpha(3:m+1)]/2/h^2;
 BB=diag(A1);
 BBBB=[];
  for k=1:m-2
     BBBB=blkdiag(BBBB,BB);
  end
BBBB=[BBBB;zeros(m-1,(m-1)*(m-2))];
BBBB=[zeros((m-1)^2,m-1) BBBB];

 C2=-1/h/2*ones(1,m-1)-[alpha(1:m-1)+alpha(2:m)]/2/h^2;
  CC=diag(C2);
   CCCC=[];
  for k=1:m-2
     CCCC=blkdiag(CCCC,CC);
  end
CCCC=[CCCC zeros((m-1)*(m-2),m-1)];
CCCC=[zeros(m-1,(m-1)^2) ; CCCC];
  
  DDDD=AAAA+BBBB+CCCC;
for k=1:n
    ff=[];
   for j=1:m-1
   ff=[ff 1/t*U(j+1,2:m)+[feval(f,h:h:1-h,h*j,k*t)]];
   end
xxx=DDDD\ff';
    xxx=reshape(xxx,m-1,m-1);
    U(2:m,2:m)=xxx';
    for j=1:m+1
   Ujing(j,1:m+1)=feval(jingqueu,(j-1)*h,x,t*k);
    end
    zuidaerr(k)=max(max(abs(U-Ujing)))
    
end
err=abs(U-Ujing);
%显示四位有效数
%zuidaERR=vpa(zuidaERR,4);
%k=n/10+1:n/10:n+1;
%U=[U(k,6)];
%U=vpa(U,7);
%ERR=[ERR(k,6)];
%ERR=vpa(ERR,4);