z=1.2;​​epsin=5.7;​​ωp=9;​​epsd=25;​​γ=0.1;​​γ0=10^-4;​​ω0=1.2;​​c=3.*10^17;​​h=6.582*10^-16;​​​​epsm[ω_]:=epsin-(ωp)^2/(ω^2+I*γ*ω);​​r[q_,ω_]:=(epsd*Sqrt[epsm[ω]/epsd-q^2]-epsm[ω]*Sqrt[1.-q^2])/(epsd*Sqrt[epsm[ω]/epsd-q^2]+epsm[ω]*Sqrt[1.-q^2]);​​​​de[ω_?NumericQ]:=Re[NIntegrate[(3*(γ0)*Sqrt[epsd]ω^3)/(4*π(ω0)^3)*(q^3*(1-r[q,ω]*Exp[2*I*(ω*Sqrt[1.-q^2]*Sqrt[epsd])/(c*h)*z]))/Sqrt[1-q^2],{q,0,100},MaxRecursion->20,Method->"LocalAdaptive"]];​​​​ak[s_?NumberQ]:=NIntegrate[de[ω]/(s+I*(ω-ω0)),{ω,0,100},MaxRecursion->30,Method->"LocalAdaptive"];​​​​InverseLaplaceTransform[1/(s+ak[s]),s,t];​​Plot[Abs[%]^2,{t,0,100}]
In[]:=
Abs[%]^2/.t->0.2
Out[]=
2
AbsInverseLaplaceTransform
1
s+ak[s]
,s,0.2