Return to Question

EDIT 2: Defining the integration of R as a function outside the NDsolve-block also seems to work.

(*DEFINE ALL THE CONSTANTS*)

Clear[ΩR, αS, κR, κS, γR, sar, c0];
ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                       (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;                       (*mechanical frequency in Hz*)
γR = 2*Pi*400;                              (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;                          (*coupling rate*)
κR = κS*((c0*ΩR)/ω0)                        (*coupling rate*)  


R[t_, P_, S_] := I κR Exp[-γR*t/2]* Integrate[P*Conjugate[S]*Exp[γR*d/2], {d, 0, t}]

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] == 
    I κS*(R[t, sP[z, t], sSt[z, t]])*sSt[z, t]}, {D[sSt[z, t], 
     z] == I κS*Conjugate[R[t, sP[z, t], sSt[z, t]]]*sP[z, t]}}
vars = {sP[z, t],sSt[z, t]};(*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};(*boundary conditions for input fields at z=0*);
sol = NDSolve[{eqs, bcs}, vars, {t, 0, 0.003}, {z, 0, 0.1}]

(* GENERATE PLOTS *)
 
DensityPlot[sP[z, t]^2 /. sol, {t, 0, 0.003}, {z, 0, 0.1}, 
 ColorFunction -> "Rainbow", 
 FrameLabel -> {"Time (s)", "Propagation distance (m)"}, 
 PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 80]

DensityPlot[sSt[z, t]^2 /. sol, {t, 0, 0.003}, {z, 0, 0.1}, 
 ColorFunction -> "Rainbow", 
 FrameLabel -> {"Time (s)", "Propagation distance (m)"}, 
 PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 80]

DensityPlot[
 Abs[R[t, sP[z, t], sSt[z, t]]] /. sol, {t, 0, 0.003}, {z, 0, 0.1}, 
 ColorFunction -> "Rainbow", 
 FrameLabel -> {"Time (s)", "Propagation distance (m)"}, 
 PlotRange -> All, PlotLegends -> Automatic, PlotPoints -> 40]

    (*DEFINE ALL THE CONSTANTS*)

Clear[\[CapitalOmega]R, \[Alpha]S, \[Kappa]R, \[Kappa]S, \[Gamma]R, 
  sar, c0];
\[Omega]0 = 
 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                   (*speed of light in vacuum*)
\[CapitalOmega]R = 2*Pi*5.611*10^6;     (*mechanical frequency in Hz*)
\[Gamma]R = 2*Pi*400;                   (*damping rate/linewidth in Hz*)
\[Kappa]S = -219.666*0.1;               (*coupling rate*)
\[Kappa]R = \[Kappa]S*((c0*\[CapitalOmega]R)/\[Omega]0)     (*coupling \
rate*)

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] == 
    I \[Kappa]S*(I \[Kappa]R Exp[-\[Gamma]R*t/2]*
       Integrate[
        sP[z, d]*Conjugate[sSt[z, d]]*Exp[\[Gamma]R*d/2], {d, 0, t}])*
     sSt[z, t]},
  {D[sSt[z, t], z] == 
    I \[Kappa]S*
     Conjugate[
      I \[Kappa]R Exp[-\[Gamma]R*t/2]*
       Integrate[
        sP[z, d]*Conjugate[sSt[z, d]]*Exp[\[Gamma]R*d/2], {d, 0, t}]]*
     sP[z, t]}}
vars = {sP[z, t], sSt[z, t]};      (*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], 
  sSt[0, t] == 
   Sqrt[0.1]};      (*boundary conditions for input fields at z=0*);
NDSolve[{eqs, bcs}, vars, {t, 0, 0.02}, {z, 0, 0.2}]

    (*DEFINE ALL THE CONSTANTS*)

Clear[ΩR, αS, κR, κS, γR, sar, c0];
ω0 = 2*Pi*2.99792458*10^8/(1.55*10^(-6));   (*optical frequency in Hz*)
c0 = 2.99792458*10^8;                       (*speed of light in vacuum*)
ΩR = 2*Pi*5.611*10^6;                       (*mechanical frequency in Hz*)
γR = 2*Pi*400;                              (*damping rate/linewidth in Hz*)
κS = -219.666*0.1;                          (*coupling rate*)
κR = κS*((c0*ΩR)/ω0)                        (*coupling rate*)  

(*SET UP PDEs TO SOLVE*)

Clear[sP, sSt, eqs, vars, bcs];
eqs = {{D[sP[z, t], z] == I κS*(I κR Exp[-γR*t/2]*
    Integrate[sP[z, d]*Conjugate[sSt[z, d]]*Exp[γR*d/2], {d, 0, t}])*sSt[z, t]},
       {D[sSt[z, t], z] == I κS*Conjugate[I κR Exp[-γR*t/2]*
    Integrate[sP[z, d]*Conjugate[sSt[z, d]]*Exp[γR*d/2], {d, 0, t}]]*sP[z, t]}}
vars = {sP[z, t], sSt[z, t]};      (*variables to solve for*)
bcs = {sP[0, t] == Sqrt[0.9], sSt[0, t] == Sqrt[0.1]};   (*boundary conditions at z=0*);
NDSolve[{eqs, bcs}, vars, {t, 0, 0.02}, {z, 0, 0.2}]