nonperturbative NLO response

clipped sine wave

amplifier response is often nonlinear
(Subscript[S, out])((Subscript[S, in]))=ArcTan[(Subscript[S, in])]

Plot[ArcTan[s],{s,-5,5}]

Clear[a];
s[t_,a_]:=a Sin[t];
Manipulate[
Plot[ArcTan[s[t,a]],{t,-10,10},PlotStyle->{Red,Thickness[lineThick]}],
{a,0.01,5}]

potentials

chi2: cubic component

a2=1;
a3=0.1;
u2[x_]:=a2 (x^2);
u3[x_]:=a3 (x^3);
Plot[{u2[x],u3[x],u2[x]+u3[x]},{x,-3,3},PlotStyle->{{Blue,Thickness[lineThick],Dashing->0.02},{Purple,Thickness[lineThick],Dashing->0.02},{Red,Thickness[lineThick]}}]

chi3: quartic component

a2=1;
a4=-0.2;
u2[x_]:=a2 (x^2);
u4[x_]:=a4 (x^4);
Plot[{u2[x],u4[x],u2[x]+u4[x]},{x,-3,3},PlotStyle->{{Blue,Thickness[lineThick],Dashing->0.02},{Purple,Thickness[lineThick],Dashing->0.02},{Red,Thickness[lineThick]}}]

adding wave components

First let's add two sine waves, at the fundamental and its second harmonic to see what the result is:

lineThick=0.004;
axesThick=0.003;

ω=2 π;
a=0.2`;
Plot[{Cos[ω t],a Cos[2 ω t],Cos[ω t]+a Cos[2 ω t]},{t,0,4},PlotStyle->{{Blue,Thickness[lineThick]},{Purple,Thickness[lineThick]},{Red,Thickness[lineThick]}}, AxesStyle->Thickness[axesThick]]

What you see is that the amplitude is asymmetric (+x is different from -x), and that there is actually a DC (constant) shift to the whole result.

Now let's see what adding the fundamental to the third harmonic looks like:

ω=2 π;
a=0.3`;
Plot[{Cos[ω t],a Cos[3 ω t],Cos[ω t]+a Cos[3 ω t]},{t,0,4},PlotStyle->{{Blue,Thickness[lineThick]},{Blue,Thickness[lineThick]},{Red,Thickness[lineThick]}}, AxesStyle->Thickness[axesThick]]

third-order nonlinearity

time domain

single frequency input, CW, step response

b = 0, linear solution

ω0=6;
ω=1;
γ=0.1;
b=0;
c=1;
tMax=300;
solution=NDSolve[{x''[t]+2 γ x'[t]+((ω0)^2) x[t]-b ((x[t])^3)==-c  Cos[ω t],x[0]==0,x'[0]==0},x,{t,0,tMax},MaxSteps->(10^6)];
xLin[t_]=x[t]/.solution[[1]];
Plot[xLin[t],{t,0,tMax}]

b > 0, include 3rd order nonlinearity

ω0=6;
ω=1;
γ=0.1;
b=0.4;
c=1;
tMax=300;
solution=NDSolve[{x''[t]+2 γ x'[t]+((ω0)^2) x[t]-b ((x[t])^3)==-c  Cos[ω t],x[0]==0,x'[0]==0},x,{t,0,tMax},MaxSteps->(10^6)];
xPert[t_]=x[t]/.solution[[1]];
Plot[xPert[t],{t,0,tMax}]

Plot[xPert[t]-xLin[t],{t,0,tMax}]

Let's zoom in to a region after the transients:

Plot[{(xPert[t]-xLin[t])/(4 (10^(-7))),xLin[t]/0.03},{t,250,tMax},PlotStyle->{{Red,Thickness[0.004]},{Blue,Thickness[0.004]}}]

pulsed input

sGauss[t_,tau_,nG_]:=Exp[-(((t/tau))^(2 nG))];

nonlinear solution
first run with ω0=5.8 (moving resonance away from any harmonic)
then with ω0=5.1 (intermediate resonance enhances higher orders)

ω0=5.1;
ω=1;
γ=0.1;
b=0.4;
c=10;
tp=200;
nSG=2;
tCtr=3 tp;
delT=(tp/10);
delT=tp;

tMax=3 tCtr;
solution=NDSolve[{x''[t]+2 γ x'[t]+((ω0)^2) x[t]-b ((x[t])^3)==-c sGauss[t-tCtr,tp,nSG] Cos[ω t],x[0]==0,x'[0]==0},x,{t,0,tMax},MaxSteps->(10^6)];
xPert[t_]=x[t]/.solution[[1]];
Plot[xPert[t],{t,tCtr-delT,tCtr+delT}]

linear solution (zero out nonlinearity)

b=0;
solution=NDSolve[{x''[t]+2 γ x'[t]+((ω0)^2) x[t]-b ((x[t])^3)==-c sGauss[t-tCtr,tp,nSG] Cos[ω t],x[0]==0,x'[0]==0},x,{t,0,tMax},MaxSteps->(10^6)];
xLin[t_]=x[t]/.solution[[1]];
Plot[xLin[t],{t,tCtr-delT,tCtr+delT}]

We can't see much difference here, so subtract off the linear solution:

Plot[xPert[t]-xLin[t],{t,tCtr-delT,tCtr+delT}]

This residual component clearly has multiple frequency components. One component is oscillating right at the fundamental. The next term is at the third harmonic.
1.00007

Frequency domain analysis 
calculate the Fourier transform of the nonlinear response: 
first sample, then use FFT

define sampling grids
numerical FT runs fastest if npts = power of 2, this makes Mathematica use the Fast-Fourier Transform algorithm
We have to use a fine sampling grid - must sample at 2x the maximum frequency we want to see in the transform. This can be substantially higher sampling rate than the fundamental driving frequency.

nptsFFT=(2^15);
Print["number of pts: ", nptsFFT];
(* dt=(2 π/ω0)/10; fixed dt, or take max range of solution and divide by npts *)
dt = tMax/nptsFFT;
tLst=Table[j dt,{j,0,nptsFFT-1}];
Print["dt: ", dt//N," time range: ",tLst[[1]]," to ",tLst[[nptsFFT]] //N];

frequency grid
Note connection between the frequency spacing and the total time range of sample

dω=2 π/(nptsFFT  dt);
ωLst=Table[j dω,{j,-nptsFFT/2,nptsFFT/2-1}];
Print["freq range: ",ωLst[[1]]//N," to ",ωLst[[nptsFFT]] //N];

sample function

xPertLst=Table[xPert[tLst[[j]]],{j,1,nptsFFT}];
ListPlot[Transpose[{tLst,xPertLst}],PlotRange->{{500,540},All}]

define FFT functions. 
See fft demo.nb for derivation of this method. Basically prefactors are there to conserve energy (satisfying Parseval's theorem).

FFT[fLst_,npts_,dt_]:=(Sqrt[((npts (dt^2))/(2 π))])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
InverseFFT[FLst_,npts_,dt_]:=(Sqrt[((2 π)/(npts (dt^2)))])RotateLeft[InverseFourier[RotateRight[FLst,npts/2]],npts/2];

transform response, plot on a Log scale
First two peaks are the fundamental and the 3rd harmonic. Can get higher-order harmonics by driving it harder. 

small hump underneath the harmonics is where the resonance frequency is

xPertLstω=FFT[xPertLst,nptsFFT,dt];

ωMin=0.`;ωMax=15;
ListPlot[Transpose[{ωLst,Abs[xPertLstω]}],Joined->True,PlotRange->{{ωMin,ωMax},All}]
ListPlot[Transpose[{ωLst,Log[Abs[xPertLstω]]}],Joined->True,PlotRange->{{ωMin,ωMax},All}]