Matlab:双曲方程

tic;
clear
clc
M=[10,20 40 80];%空间步数
N=2*M;%时间步数
for k=1:length(M)
h=1/M(k);%空间步长
tau=1/N(k);%时间步长
s=tau/h;%步长比
x=0:h:1;
t=0:tau:1;
y=inline('exp(x+t)','x','t');%真解函数
for i=1:length(x)
    for j=1:length(t)
        exact(i,j)=y(x(i),t(j));%真解
    end
end
u=zeros(M(k)+1,N(k)+1);%数值解内存单元
for i=2:M(k)
    u(i,1)=exp(x(i));%初值u(i,0)
    u(i,2)=exp(x(i))+tau*exp(x(i))+(tau^2/2)*exp(x(i));%第二层值u(i,1)
end
u(1,:)=exp(t);%边值u(0,t)
u(M(k)+1,:)=exp(1+t);%边值u(1,t)
phi=zeros(M(k)-1,N(k));
for i=1:N(k)
    phi(1,i)=exp(t(i));
    phi(M(k)-1,i)=exp(1+t(i));
end

A=diag((2*(1-s^2))*ones(M(k)-1,1))+diag(s^2*ones(M(k)-2,1),1)+diag(s^2*ones(M(k)-2,1),-1);
for j=2:N(k)
    u(2:M(k),j+1)=A*u(2:M(k),j)-u(2:M(k),j-1)+s^2*phi(:,j);%数值解
end
error=abs(u(2:M(k),2:end)-exact(2:M(k),2:end));%误差
error_inf(k)=max(max(error));
%error=exact-u;
subplot(1,2,1)
[X,Y]=meshgrid(t(end:-1:2),x(2:M(k)));
mesh(X,Y,error);
xlabel('t');
ylabel('x');
zlabel('error');
grid on 
% pan on;
% zoom on;
pause(0.05)
hold on
end
legend('h=1/10,tau=1/20','h=1/20,tau=1/40','h=1/40,tau=1/80','h=1/80,tau=1/160');
title('Numerical Result')
for p=1:length(M)-1
H=error_inf(p)/error_inf(p+1);
NORM(p)=log2(H);
end
subplot(1,2,2)
plot(1:length(M)-1,NORM,'-bp')
ylabel('误差阶数');
text(2,2,'这就是误差阶数!!!!!')
grid on
axis square
toc; 

效果图：