In this appendix, we show how to use the model based on Schröder et al. to simulate the DFWM mode-locked fiber laser. The following code is based on [65] and written in Mathematica 9.0:
 (* This code solves the NLS equation with the split-step Fourier method
based on Govind P. Agrawal in March 2005 for the NLFO book
*)

ClearAll["Global`*"]

(*---Specify input parameters*)
distance = 150.;(*Enter fiber length (in units of L_c)=*)
kappa = -0.001 ;(*Normalized 2nd-dispersion: kappa=beta2*f^2*L/2): \
+ve for normal,-ve for anomalous*)
sigma = 0.; (*Normalized 3rd-dispersion: sigma=beta3*f^3*L/6 *)
G = 1.0;(*small signal gain of Ramam Amp: G=g*L *)
Is = 1.; (*gain saturation parameter*)
alpha = 0.4; (*Normalized fiber amplitude absorption coeff: alpha=l*L*)


n = 2.^0.5;(*Nonlinear parameter n=') \
sqrt(L_D/L_NL)=sqrt(gamma*P0*T0^2/|beta2|) or QT: n=kappa^0.5*)

(*---Specify filter parameters*)
bdwidth = 2. Pi*6.;
delta = 2. Pi*0.5;
a = Log[Sqrt[0.95]];
perta = 0.3;
pertfsr = 0.17;
T = 0.2;
t = Sqrt[T];
r = I*Sqrt[1 - T];

(*---Specify input parameters*)
mshape = -1.;(*m=0 for sech,m>0 for super-Gaussian=*)
chirp0 = 0. ;(*% input pulse chirp (default value)*)

(*
P = 1/(gamma*L); (P is the ref peak power);  
uu = A/sqrt(P); 
z = z0/L_c; (z0 is the real length, L_c is the cavity length);
tau = f*t; (t is the time of reference traveling frame);
*)

(*---set simulation parameters*)
nt = 2^13; (*% FFT points  (powers of 2)*)
Tmax = 100.; (*(half) window size*)
stepno = 1*Round[20*distance*n^2];(*No.of z steps to*)
dz = distance/stepno;(*step size in z*)
dtau = (2.*Tmax)/nt ;(*step size in tau*)

Twin = 5.;
fmax = (1./(2. Tmax))*nt/2.;
fwin = 5.;

filterz = 0.5;
plotz = 5;

(*---tau and omega arrays*)
tau = N[Range[-nt/2, nt/2 - 1]]*dtau ;(*temporal grid*)
omega = (1. \[Pi]/Tmax)*
   N[Join[Range[0, nt/2 - 1], 
     Range[-nt/2, -1]]];(*[(0:nt/2-1) (-nt/2:-1)]*)
(*frequency grid*)
delaytau = dtau*Range[-Round[Twin/dtau], Round[Twin/dtau]];



(*Input Field profile*)
If[mshape == 0,
  (*;% soliton*)
  uu = Sech[tau]*Exp[-0.5 I*chirp0*tau^2.],
  If[ mshape > 0 ,
   (* super-Gaussian*)
   uu = Exp[-0.5*(1. + 1. I*chirp0)*tau^(2.*mshape)],
   (*White noise*)
   uu = (RandomReal[NormalDistribution[0, 1], nt] + 
       I*RandomReal[NormalDistribution[0, 1], nt])*Sqrt[1./2.]
   ]
  ];

temp = RotateRight[
   InverseFourier[uu, FourierParameters -> {1, -1}]*(nt*dtau)/
     Sqrt[2.*Pi], nt/2.];
tempomega = RotateRight[omega, nt/2];


(*---store dispersive phase shifts to speedup code*)
dispersion = 
  Exp[(-alpha + I*kappa*omega^2. + I*sigma*omega^3.)*
    dz];(*% nonlinear phase factor*)
(*comb filter type*)
(*original comb filter + BPF*)
(*filtert= \
Exp[-omega^2./bdwidth^2.]*(t^2)/(1-r^2*Exp[-I*(omega+delta)+a]);*)
(*perturbated comb filter + BPF*)
filtert = Exp[-omega^2./bdwidth^2. - perta*Sin[0.5
         *omega/pertfsr]^2]*(t^2)/(1 - 
      r^2*Exp[-I*(omega + delta) + a]);

p6 = ListLinePlot[
   Transpose[{tempomega/(2. Pi), 
     RotateRight[10.*Log10[Abs[filtert]^2], nt/2]}], 
   PlotRange -> {{-fwin, fwin}, {-30, 10}}, Frame -> True, 
   ImageSize -> Automatic];

(*%*********[Beginning of MAIN Loop]***********
% scheme:1/2N\[Rule]D\[Rule]1/2N;first half step nonlinear*)
temp = uu*
   Exp[(I*Abs[uu]^2. + G/(1.0 + Abs[uu]^2./Is))*dz/2.];(*% note hhz/2*)

starttime = SessionTime[];
timeused = SessionTime[] - starttime;
z = 0;

Monitor[ (* Realtime monitoring the simulation progress*)
 For[ i = 1, i <= stepno, i++,
  If[Round[(Mod[z, 1] - filterz )/dz] == 0 ,
    ftemp = 
     InverseFourier[temp, FourierParameters -> {1, -1}]*filtert*
      dispersion;
    ftemp = 
     InverseFourier[temp, FourierParameters -> {1, -1}]*dispersion;
    ]
   uu = Fourier[ftemp, FourierParameters -> {1, -1}];
  temp = uu*Exp[(I*Abs[uu]^2. + G/(1.0 + Abs[uu]^2./Is))*dz];
  z = z + dz;
  If[Round[Mod[z, plotz ]/dz] == 0 || 
    Round[(Mod[z, plotz ] - plotz )/dz] == 0 ,
   timeused = SessionTime[] - starttime;
   p1 = ListLinePlot[Transpose[{tau, Abs[temp]^2.}], 
     PlotRange -> {{-Twin, Twin}, All}, Frame -> True, 
     FrameLabel -> {{Null, Null}, {Null, 
        StringJoin[{"i = ", ToString[i], ", z = ", ToString[z]}]}}
     ];
   ftemp0 = RotateRight[ftemp (nt*dtau)/Sqrt[2*Pi], nt/2];
   p2 = ListLinePlot[
     Transpose[{tempomega/(2. Pi), 10.*Log10[Abs[ftemp0]^2.]}], 
     PlotRange -> {{-fwin, fwin}, All}, Frame -> True, 
     FrameLabel -> {{Null, Null}, {Null, 
        StringJoin[{" Time = ", ToString[timeused]}]}}
     ];
   p3 = ListLinePlot[Transpose[{tau, Abs[temp]^2}], 
     PlotRange -> {{-Tmax, Tmax}, All}, Frame -> True, 
     FrameLabel -> {{Null, Null}, {Null, 
        StringJoin[{"i = ", ToString[i], ", z = ", ToString[z]}]}}
     ];
   p4 = ListLinePlot[
     Transpose[{tempomega/(2. Pi), 10.*Log10[Abs[ftemp0]^2.]}], 
     PlotRange -> {{-fmax, fmax}, All}, Frame -> True, 
     FrameLabel -> {{Null, Null}, {Null, 
        StringJoin[{" Time = ", ToString[timeused]}]}}
     ];
   autocorr0 = 
    RotateRight[
     InverseFourier[
      Fourier[Abs[temp]^2]*Conjugate[Fourier[Abs[temp]^2]]], nt/2];
   p5 = ListLinePlot[Transpose[{tau, autocorr0/Max[autocorr0]}], 
     PlotRange -> {{-Twin, Twin}, All}, Frame -> True];
   p = GraphicsGrid[{{p1, p2}, {p3, p4}, {p5, p6}}, ImageSize -> Full]
   ];
  ],
 p
 ]

uu = temp*
   Exp[(I*Abs[uu]^2. + G/(1.0 + Abs[uu]^2./Is))*
     dz/2.];(*% Final field*)
temp = RotateRight[
   InverseFourier[uu, FourierParameters -> {1, -1}]*(nt*dtau)/
     Sqrt[2*Pi], nt/2];

(*%Final spectrum*)
(*%***************[End of MAIN Loop]***************)
p1 = ListLinePlot[Transpose[{tau, Abs[uu]^2}], 
   PlotRange -> {{-Twin, Twin}, All}, Frame -> True, 
   FrameLabel -> {{Null, Null}, {Null, 
      StringJoin[{"i = ", ToString[i], ", z = ", ToString[z]}]}}
   ];
p2 = ListLinePlot[
   Transpose[{tempomega/(2. Pi), 10.*Log10[Abs[temp]^2.]}], 
   PlotRange -> {{-fwin, fwin}, All}, Frame -> True, 
   FrameLabel -> {{Null, Null}, {Null, 
      StringJoin[{" Time = ", ToString[timeused]}]}}
   ];
p3 = ListLinePlot[Transpose[{tau, Abs[uu]^2}], 
   PlotRange -> {{-Tmax, Tmax}, All}, Frame -> True, 
   FrameLabel -> {{Null, Null}, {Null, 
      StringJoin[{"i = ", ToString[i], ", z = ", ToString[z]}]}}
   ];
p4 = ListLinePlot[
   Transpose[{tempomega/(2. Pi), 10.*Log10[Abs[temp]^2.]}], 
   PlotRange -> {{-fmax, fmax}, All}, Frame -> True, 
   FrameLabel -> {{Null, Null}, {Null, 
      StringJoin[{" Time = ", ToString[timeused]}]}}
   ];
autocorr0 = 
  RotateRight[
   InverseFourier[Fourier[Abs[uu]^2]*Conjugate[Fourier[Abs[uu]^2]]], 
   nt/2];
p5 = ListLinePlot[Transpose[{tau, autocorr0/Max[autocorr0]}], 
   PlotRange -> {{-Twin, Twin}, All}, Frame -> True];
p = GraphicsGrid[{{p1, p2}, {p3, p4}, {p5, p6}}, ImageSize -> Full]

(*Exporting results*)
fname = "2nd_01(0.18)"; (*file name*)
Export[StringJoin[{fname , "_time.dat"}], 
  Transpose[{tau, Re[uu], Im[uu], Abs[uu]^2}]];
Export[StringJoin[{fname , "_freq.dat"}], 
  Transpose[{tempomega/(2. Pi), Re[temp], Im[temp], Abs[temp]^2, 
    RotateRight[Abs[filtert]^2, nt/2]}]];
Export[StringJoin[{fname , "_autocorr.dat"}], 
  Transpose[{tau, autocorr0, autocorr0/Max[autocorr0]}]];