clear;
Lx=1; Ly=1; T=0.05;
Nx=31; Nt=200;
Dx=Lx/(Nx-1); Dt=T/Nt;
Tn=zeros(Nx,Nx);
x=linspace(0,Lx,Nx);
y=linspace(0,Lx,Nx);
[X,Y]=meshgrid(x,y);
%начальные условия
for n=1:Nx
for m=1:Nx
Tn(n,m)=0.1*sin(pi*Dx*(n-1))*sin(pi*Dx*(m-1))*exp(i*(Dx*(n-1))^2); %0.1*sin(pi*x)*sin(p*y)*exp(i*x^2)
end
end
%продублировал граничные условия на всякий, если начальные условия будут плохие
Tn(1,:)=0;
Tn(Nx,:)=0;
Tn(:,1)=0;
Tn(:,Nx)=0;
%основной цикл
for n=1:Nt
Tc=Tn;
for q=2:Nx-1
for m=2:Nx-1
Tn(q,m)=((1+i/Dt)*Tc(q,m))/(i/Dt)+((Tc(q,m+1)+Tc(q+1,m)-4*Tc(q,m)+Tc(q,m-1)+Tc(q-1,m))/(Dx^2))/(i/Dt);
end
end
%граничные штуки
Tn(1,:)=0;
Tn(Nx,:)=0;
Tn(:,1)=0;
Tn(:,Nx)=0;
figure(1)
mesh(x,y,abs(Tn));
end