(*Response of a SH oscillator to a periodic half-wave rectified sine input*)
ClearAll["@"]
h[t_]=UnitStep[t];

(*Mathematica can solve if we truncate the input to a finite sum*)
(*E.g. summing to n=10 gives the exact solution for t < 11Pi*)
de = x''[t]+4x[t] == Sum[Sin[t-n*Pi]*h[t-n*Pi], {n,0,10}];
ic1 = x[0]==0;
ic2 = x'[0]==0;
soln = FullSimplify[DSolve[{de,ic1,ic2}, x[t], t]]
Plot[x[t] /. soln, {t,0,33}]
