这两个程序用于仿真 图4和图5.这里的方程用新的归一化方程，不是孤子归一化方程。也就是论文中的方程3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%图4%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  study input peak power vursus output total power for different Lcd, P~P CURVES 
%  Simulate Fig.4 Using wangs' equations;

lcd=0.1;
M1=43,
N =512;                         % Number of Fourier modes (Time domain sampling points)
dz =0.5*0.0031415926/4;         % space step, make sure nonlinear<0.05
M3=4000;
T =40;                          % length of time:T*T0, it can.t be too big ot too small, it affect accuracy
T0=60;                          % input pulse width
dt = T/N;                       % time step
n = [-N/2:1:N/2-1]';            % Indices
t = n.*dt; 
ww = 4*n.*n*pi*pi/T/T;          % Square of frequency. Note i^2=-1.
w=2*pi*n./T;
www=w.*ww;
g1=-i*ww./(2*lcd);
g2=-i*ww./(2*lcd);         %w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0;
P1=0;
P2=1;
P=0;
for m1=1:1:M1
p=0.5+0.03*m1,
%s10=p.*sech(p.*t).*exp(0.5*sqrt(-1.)*t.*t);
s10=p.*sech(t);                                        % input
%s10=p.*exp(-0.5*(1+sqrt(-1.))*t.*t);
%s10=p.*exp(-t.*t./(2*1.146*1.146));                  
%s10=p.*exp(-t.*t./2);
s1=s10;
s20=0*s10;
s2=s20;
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1)))); %energy in waveguide 1
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1)))); %energy in waveguide 2
for m3 = 1:1:M3                                    % Start space evolution
   s1 = exp(dz*i*4*(abs(s1).*abs(s1))).*s1;        % Solve nonlinear part of NLS
   s2 = exp(dz*i*4*(abs(s2).*abs(s2))).*s2;
   sca1 = fftshift(fft(s1));                       % Take Fourier transform
   sca2 = fftshift(fft(s2));
   sc2=exp(g2.*dz).*(sca2+i*sca1.*dz);
   sc1=exp(g1.*dz).*(sca1+i*sca2.*dz);  % frequency domain phase shift   % Advance in Fourier space
   s2 = ifft(fftshift(sc2));                       % Return to physical space
   s1 = ifft(fftshift(sc1));
 end
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1))));
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1))));
   P1=[P1 p1/p10];
   P2=[P2 p2/p10];
   P=[P p*p];
end
figure(1)
plot(P,P1, P,P2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%Fig5%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  This Matlab script file solves the nonlinear Schrodinger equation
%  simulate fig.5 using wangs' equation                     
N = 1024;                                   % Time domain sampling points
lcd=0.2;                                    % product of dispersion length and coupling coefficient: CL_D
dz =0.00314/8;                              % space step, make sure nonlinear<0.05
M3=1000;                                    % M3*dz is the length of coupler
J =200;                                     % Steps between output of space
T =50;                                      % length of time:T*T0
MN3=0;
dt = T/N;                                    % time step
n = [-N/2:1:N/2-1]';                         % Indices
t = n.*dt; 
ww = 4*n.*n*pi*pi/T/T;                       % Square of frequency. Note i^2=-1.
w=2*pi*n./T;
g1=-i*ww./(2*lcd);
g2=-i*ww./(2*lcd);             %w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./T
s1=sech(t);                                  % input
s2=s1.*0.;
s10=s1;
s20=s10*0.;
S1=s1;
S2=s2;
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s10(1,1)*s10(1,1)))); %energy in waveguide 1
p20=dt*(sum(abs(s20').*abs(s20'))-0.5*(abs(s20(N,1)*s20(N,1))+abs(s20(1,1)*s20(1,1)))); %energy in waveguide 2
P1=p10/p10;
P2=p20/p10;
for m3 = 1:1:M3                                  % Start space evolution
   s1 = exp(dz*i*4*(abs(s1).*abs(s1))).*s1;      % Solve nonlinear part of NLS
   s2 = exp(dz*i*4*(abs(s2).*abs(s2))).*s2;
   sca1 = fftshift(fft(s1));                     % Take Fourier transform
   sca2 = fftshift(fft(s2));
   sc2=exp(g2.*dz).*(sca2+i*sca1.*dz);
   sc1=exp(g1.*dz).*(sca1+i*sca2.*dz);           % frequency domain phase shift 
                                                 % Advance in Fourier space
   s2 = ifft(fftshift(sc2));                     % Return to physical space
   s1 = ifft(fftshift(sc1));
   p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1))));
   p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1))));
if rem(m3,J) == 0                                % Save output every J steps.
   S1 = [S1 s1];                                 % put solutions in U array
   S2=[S2 s2];
   P1=[P1,p1/p10];
   P2=[P2,p2/p10];
   MN3=[MN3 m3];
   z3=dz*MN3';                                   % output location
 end
end
%figure(3)
%waterfall(t',z3',abs(S1').*abs(S1'))  %t' is 1xn, z' is 1xm, and U1' is mxn
%figure(5)
%waterfall(t',z3',abs(S2').*abs(S2'))  %t' is 1xn, z' is 1xm, and U1' is mxn
figure(1)
plot(t, abs(s2').*abs(s2'),'--', t, abs(s1').*abs(s1'))