高斯光束的非线性薛定谔方程

% idu/dz = - 1/2 * ( d^2u/dx^2 + d^2u/dy^2 - beta d^2u/dt^2) - [1-sg*R(x,y]*|u|^2*u - p*R(x,y)*u
clear all;clc;

% 计算厄米特多项式
syms H1 H2 Hn hr
n=1;%input('请输入厄米特多项式阶数:');
Hn(1:n+1)=nan;
if n==0
Hn(1)=1;
else if n==1
Hn(2)=2*hr;
else
Hn(1)=1;Hn(2)=2*hr;
for ni=3:n+1
H1=Hn(ni-2);
H2=Hn(ni-1);
Hn(ni)=2*hr*H2-2*(ni-2)*H1;
end
end
end
H=expand(Hn(n+1));
% END

%---初始化
z = 4;
np=100; Xmax=8; Ymax=8; Tmax=6;
step_num = z*10;
deltaz = z/step_num;
Z=0:deltaz:z-deltaz;
dd = (2*Xmax)/np;
dt = (2*Tmax)/np;
x = (-np/2:np/2-1)*dd; y = x;
t = (-np/2:np/2-1)*dt;
kx = (pi/Xmax)*[(0:np/2-1) (-np/2:-1)]; ky = kx;
kt = (pi/Tmax)*[(0:np/2-1) (-np/2:-1)];
[X Y T]=meshgrid(x,y,t);
[KX KY KT]=meshgrid(kx,ky,kt);

uu = 0.05*exp(-(X.^2+Y.^2+T.^2)/3);

%---厄米特多项式赋值
Jx=subs(H,X);
Jy=subs(H,Y);
V=Jx.*Jy.*exp(-(X.^2+Y.^2)./4);
V=V/max(max(max(V)));
%{
figure
surf(x,y,abs(V(:,:,4)));
shading interp
%}
%---end

beta = 12;
sg = 0.2;
p = 0.5;
dispersion =exp( -1i/2*(KX.^2+KY.^2-beta*KT.^2)*deltaz);


figure
surf(x,y,abs(uu(:,:,50)).^2);
shading interp
xlabel('X');
ylabel('Y');
zlabel('|u|^2');
title('Input Pulse');

%uz(1:np,1:np,1:z) = nan;
for n=1:step_num
display(sprintf('%d',n));
hhz = exp(-1i*((1-sg.*V).*abs(uu).^2+p*V)*deltaz);
A1 = uu.*hhz;
A2 = fftn(A1);
A3 = A2.*dispersion;
A4 = ifftn(A3);
uu = A4;
%{
if ~mod(step_num,z) && (step_num>9)
uz(:,:,n)=uu(:,:,50);
end
%}
end

figure
surf(x,y,abs(uu(:,:,50)).^2);
shading interp
xlabel('X');
ylabel('Y');
zlabel('|u|^2');
title('Output Pulse');


figure
isosurface(x,y,t,uu,0.008);
axis equal;
colormap([0 0 1]);