function drying
% 初始化参数
clear;
d1=2;d2=2;
r=1;k=500000;mmu1=0.1;alpha=1;mmu2=0.3;A=0.2;
tau=0.4;
l=30;%区间长度为l*pi

N=20;%N是tau/h(h为时间上划分的步长)
h=tau/N;
K=6000;%时间方向计算步数
T=h*K;%时间最值为T
p=K+N+1;%时间点总数
M=100;%Delta{x}为空间上的划分,M=pi/Delta{x}
dx=l/M;

%初始化，求u，v在每个网格的节点上的值
u=zeros(M+1,p);
v=zeros(M+1,p);

%赋初值，
for i=1:M+1
    for j=1:N+1
 %        u(i,j)=0.2+0.001*sin((i-1)*dx);
 %        v(i,j)=0.02+0.02*sin((i-1)*dx);
        u(i,j)=0.5;
        v(i,j)=4.5;
%         u(i,j)=beta*sin((i-1)*dx);
%         v(i,j)=(a-b*beta)/c*sin((i-1)*dx);
    end 
end

% 求下一组节点的值
f=zeros(M+1,p);
g=zeros(M+1,p);
for j=N+1:p-1 
    for i=2:M
        f(i,j)=r*u(i,j)-v(i,j-N)/(1+alpha*v(i,j-N))-mmu1*u(i,j);
        g(i,j)=beta*u(i,j)*v(i,j-N)-(mmu2-A)*v(i,j);
        u(i,j+1)=u(i,j)+h*d1*(u(i+1,j)-2*u(i,j)+u(i-1,j))/dx/dx+h*f(i,j);
        v(i,j+1)=v(i,j)+h*d2*(v(i+1,j)-2*v(i,j)+v(i-1,j))/dx/dx+h*g(i,j);
    end
   % 赋边值，x=0，l*pi时，即x方向在第1,M+1的值neuman
    u(1,j+1)=u(2,j+1);  
    u(M+1,j+1)=u(M,j+1);
    v(1,j+1)=v(2,j+1);  
    v(M+1,j+1)=v(M,j+1);
end

%作图
z1=zeros(p,M+1);
z2=zeros(p,M+1);
for i=1:M+1
    for j=1:p
        z1(j,i)=u(i,j);
        z2(j,i)=v(i,j);
    end
end

x=0:dx:l;
t=-r:h:T;
[xx,tt]=meshgrid(x,t);
figure
% surf(x,t,z1);
mesh(xx,tt,z1);
figure
mesh(xx,tt,z2+0.3);
% surf(x,t,z1);

有所更改，自取⑧