Clear[f,g,h,t,tp,r,ω0,ω,Kernel,b1,Ψ1,Ψ2,β,ψ,MAX];

f[t_]=-I 2β Cos[ω*t]Exp[-I ω0*t];
g[tp_,t_] = DiracDelta[tp-t];
h[t_]=-I 2β Cos[ω*t]Exp[I*ω0*t];

ω0=r ω; (*r = 1 on resonnance. Leave r unevaluated for off-resonance work. Put r=1 NOW if working on resonance as the limit r->1 in the kernel is non-trivial and taking it after the computations leads to erroneous 1/0*)

r=1;(*Working on resonance for the example, we put r=1. Remove this to work off resonance*)

(*Calculation of the kernel K_up  in the interaction picture. Passing to the interaction picture is not required*)
(*All parameters are left as formal variables so as to get the universally valid analytical expansion of the path-sum solution. 
You can speed up all calculations by giving values to the parameters now. Use exact values to keep analytical calculations, e.g. write b = (2+3/10) ω instead of 2.3 ω
	Use numerical values e.g. b=2.3 ω for further speed up
*)
Kernel[tp_,t_]= Integrate[Integrate[f[tp]*g[τ2,τ1]*h[τ1],{τ2,τ1,tp},Assumptions->{tp>0,t>=0,tp≥t}],{τ1,t,tp},Assumptions->{tp>0,t>=0,tp>t}]//FullSimplify(*This is K_up. For K_down, just exchange h and f in the formula!*)

(*Computation of G_up up to order MAX*)
Clear[Res,MAX,jj];
Res[t_]=Kernel[t,0];
MAX=4; (*Enter the maximum desired order here*)

(*This will display the number of the order reached once the calculation of this order is finished*)
Clear[jj,Store];
For[jj=1,jj<MAX+1,
Store=Integrate[Kernel[t,τ]*Res[τ],{τ,0,t},Assumptions->{t>0}];
ToExpression["Res"<>ToString[jj+1]<>"[t_]=Store"];
Clear[Res];
Res[t_]=Store;
Print[jj];
jj++;
];Clear[jj];

Res1[t_]=Kernel[t,0];

(*Iterative integration of all orders of G_up*)
For[jj=1,jj<MAX+1,
ToExpression["Isq"<>ToString[jj]<>"[t_]=Integrate[Res"<>ToString[jj]<>"[τ],{τ,0,t}]"];
Print[jj];
jj++;
];Clear[jj];

(*last iterative integration as required by the Neumann expansion of the path-sum solution*)
For[jj=1,jj<MAX+1,
ToExpression["Isq2"<>ToString[jj]<>"[t_]=Integrate[h[τ]*Isq"<>ToString[jj]<>"[τ],{τ,0,t}]"];
Print[jj];
jj++;
];Clear[jj];

(*Now we get access to the orders generated for U_{up->down}. Here I wrote only two but all are accessible up to MAX. Replace h[τ] by f[τ] if calculating U_{down->up}*)
H[t_]=Integrate[h[τ],{τ,0,t}];(*This is order 0*)

Clear[Ampl2,Ampl3,Ampl4];
Ampl2[t_]=H[t]+Sum[ToExpression["Isq2"<>ToString[jj]<>"[t]"],{jj,1,2}];(*Order 2*)
P2[t_]=Abs[Ampl2[t]]^2;(*Transition probability order 2*)

Ampl3[t_]=H[t]+Sum[ToExpression["Isq2"<>ToString[jj]<>"[t]"],{jj,1,3}];(*Order 3*)
P3[t_]=Abs[Ampl3[t]]^2;(*Transition probability order 3*)

Ampl4[t_]=H[t]+Sum[ToExpression["Isq2"<>ToString[jj]<>"[t]"],{jj,1,4}];(*Order 4*)
P4[t_]=Abs[Ampl4[t]]^2;(*Transition probability order 4*)

(*Let's look at the formulas in details, EXEMPLE SHOWN FOR P3 but works for all orders*)

Simplify[P3[t],Assumptions->{Element[{β,ω,r,t},Reals]}](*Mathematica has trouble evaluating absolute values analytically.*)

(*Instead we force it to give us the formulas*)

Ampl3c[t_]=Ampl3[t]//.{I->(-I)};(*Conjugate of amplitude*)

(* And P3 = Ampl3 * Ampl3c *)

(*Now we can get access to the fully explicit form of P, here of P2:*)
(*Function S0 is always 0*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[2ω t]](*Function S2*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[4ω t]](*Function S4*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[6ω t]](*Function S6*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[8ω t]](*Function S8*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[10ω t]](*Function S10*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[12ω t]](*Function S12*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[14ω t]](*Function S14*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Sin[16ω t]](*Function S16*)

(ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]]//.{Sin[x_]->0,Cos[y_]->0})(*For C0, we artificially put all sines and cosines to 0 since it is the part of the solution with no trigonometric factor*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[2ω t]](*Function C2*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[4ω t]](*Function C4*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[6ω t]](*Function C6*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[8ω t]](*Function C8*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[10ω t]](*Function C10*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[12ω t]](*Function C12*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[14ω t]](*Function C14*)
Coefficient[ExpToTrig[Expand[Ampl3[t]Ampl3c[t]]],Cos[16ω t]](*Function C16*)

(*Let's plot the results and compare with a fully numerical solution*)
Clear[β,Ψ1,Ψ2,t,Ea,Eb,ω,ω2,r]

(*Now we can put parameter values*)
r=1; (*If off-resonance, put another value of r than 1 here*)
ω0=1;ω=ω0*1/r;

β=0.5*ω;
TMAX=9;(*In units of ωt*)

(*Fully numerical solution for comparison. Used here is the original Hamiltonian*)
sol=NDSolve[{-I *(ω0/2)*Ψ1[t]-2 I β Cos[ω t]Ψ2[t]==D[Ψ1[t],t],-2 I β Cos[ω t]Ψ1[t]-I*(-ω0/2)*Ψ2[t]==D[Ψ2[t],t],Ψ1[0]==1,Ψ2[0]==0},{Ψ1[t],Ψ2[t]},{t,0,4TMAX}];
func[t_]=Ψ2[t]/.sol;

(*Example with orders 3 & 4 here:*)
Show[{
Plot[{Chop[P3[Wt/ω]],Chop[P4[Wt/ω]],Abs[func[Wt/ω]]^2},{Wt,0,TMAX}(*WT is in units of ω t*)
,PlotRange->All,
PlotStyle->{Directive[Dotted,Red],Directive[Blue],Directive[Dashed,Black]},Frame->True,Axes->None,FrameLabel->{"ω t","(Subscript[P, ↑↦↓])(t)"},ImageSize->600,FrameStyle->Directive[Black,24]]
},PlotRange->{0,1},ImageSize->600]
Clear[β,Ψ1,Ψ2,t,Ea,Eb,ω,ω2,r];