如下Matlab 程序计算 论文Youfa Wang and Wenfeng Wang,  Study of ultrafast pulse coupling dynamics considering retarded
nonlinear response and self-steepening effects, Joutnal of lightwave Technology, Vol.24, No.2,
pp1041-1047,中的 图2。 做一点简单改动，文中其他图，也可以计算。

% nonlinear equation including nonlinear delay, self-steepening, and higher order dispersion
%  Youfa Wang and Wenfeng Wang,  Study of ultrafast pulse coupling dynamics considering retarded
%  nonlinear response and self-steepening effects, Joutnal of lightwave Technology, Vol.24, No.2,
%  pp1041-1047,
MN3=0;
lcd=0.5;                    %lcd>0 corresponding to beita2<0, lcd<0 corresponds
%lcd=1000;
ldd=-200*lcd;
lmd=-0.00;
fr=0.18;
st=0.028;                         %st=1/(w0*t0)
%st=0;
J=400;
N =1024;                          % Number of Fourier modes (Time domain sampling points)
M3=6000;
dz =1.570796326/M3;               % space step, make sure nonlinear<0.05
T =20;                            % normalized time, t=T*T0, it can't be too big ot too small, it affect accuracy
delt=-0.0;                        % normal mismatch=(b1-b2)/2c
T0=28.4;
dt = T/N;                         % time step
n = [-N/2:1:N/2-1]';              % Indices
t = n.*dt;
%----Ranman response in fiber---------------------------
n1=[0:1:N/2]';         % the response is 0, if t1<0. it is not easy to use t. the result in funcation no difference
t1=n1*dt;
ht=(32*32+12.2*12.2)/(32*32*12.2)*T0*exp(-t1*T0/32).*sin(t1*T0/12.2);
%----------------display convolution-------------
y=[-N/2:1:N-1]';
t2=dt*y;
%--------------------------------------------------------------
ww = 4*n.*n*pi*pi/T/T;            % Square of frequency. Note i^=-1.
w=2*pi*n./T;
www=w.*ww;
g1=i*(delt/T0)*w-i*ww./(2*lcd)-i*www./(6*ldd);
g2=i*(delt/T0)*w-i*ww./(2*lcd)-i*www./(6*ldd);             %w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0;
p=1.2;
s10=p.*sech(t);
s1=s10;
s20=0*s10;
s2=s20;
S1=s1.*0; S2=s1.*0;SC1=s1.*0; SC2=s1.*0;
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
for m3 = 1:1:M3                                    % Start space evolution
   s1squre=abs(s1.*s1);
   s1cubic=s1squre.*s1;
   s2squre=abs(s2.*s2);
   s2cubic=s2squre.*s2;
   chta=conv(ht,s1squre);
   cht1=dt*chta(1:N);
   chtb=conv(ht,s2squre);
   cht2=dt*chtb(1:N);
   s1(1)=s1(1)+4*(1-fr)*i*dz*s1cubic(1)-4*(1-fr)*st*dz/dt*s1cubic(1)+4*fr*i*dz*s1(1).*cht1(1)-...
   4*fr*st/dz*s1(1).*cht1(1);
   s1(2:N)=s1(2:N)+4*(1-fr)*i*dz*s1cubic(2:N)-4*(1-fr)*st*dz/dt*(s1cubic(2:N)-...
   s1cubic(1:N-1))+4*fr*i*dz*s1(2:N).*cht1(2:N)-4*fr*st*dz/dt*(s1(2:N).*cht1(2:N)-s1(1:N-1).*cht1(1:N-1));
   s2(1)=s2(1)+4*(1-fr)*i*dz*s2cubic(1)-4*(1-fr)*st*dz/dt*s2cubic(1)+4*fr*i*dz*s2(1).*cht2(1)-...
   4*fr*st/dz*s2(1).*cht2(1);
   s2(2:N)=s2(2:N)+4*(1-fr)*i*dz*s2cubic(2:N)-4*(1-fr)*st*dz/dt*(s2cubic(2:N)-...
   s2cubic(1:N-1))+4*fr*i*dz*s2(2:N).*cht2(2:N)-4*fr*st*dz/dt*(s2(2:N).*cht2(2:N)-s2(1:N-1).*cht2(1:N-1));
   sca1 = fftshift(fft(s1));                       % Take Fourier transform
   sca2 = fftshift(fft(s2));
   sc2=exp(g2.*dz).*(sca2+i*(1+lmd*w).*sca1.*dz);
   sc1=exp(g1.*dz).*(sca1+i*(1+lmd*w).*sca2.*dz);  % frequency domain phase shift  
   s2 = ifft(fftshift(sc2));                       % Return to physical space
   s1 = ifft(fftshift(sc1));
  % if rem(m3,J) == 0                              % Save output every J steps.
  % S1 = [S1 s1];                                  % put solutions in U array
  % S2=[S2 s2];
  % SC1 = [SC1 sc1];                               % put solutions in U array
  % SC2=[SC2 sc2];
  % MN3=[MN3 m3];
  % z3=dz*MN3';                                    % output location
%end
m3,a=max(abs(cht1)),
 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+p2-p10,
%figure(3)
%waterfall(w',z3',abs(SC1').*abs(SC1'))  %t' is 1xn, z' is 1xm, and U1' is mxn
%figure(4)
%waterfall(w',z3',abs(SC2').*abs(SC2'))  %t' is 1xn, z' is 1xm, and U1' is mxn
figure(1)
plot(t,abs(s1'.*s1'), t,abs(s2'.*s2'),'.');
figure(2)
plot(t,abs(s2'.*s2'),'.');

%figure(2)
%plot(w,abs(sc1'.*sc1'), w,abs(sc2'.*sc2'));
%-----------display convolution---
%plot(t2,chta)
%------------------------------------

http://blog.163.com/opto_wang/