(************************************************)
(* Regular perturbation: y'' = Sin(x) y;
   y(0) = 1;  y'(0) = 1                        *)
(************************************************)

Remove["Global`*"]
s0 = DSolve[{p''[x]==0,p[0]==1,p'[0]==1},p[x],x][[1]];
y0 = Evaluate[p[x]/.s0]

s1 = DSolve[{q''[x]==Sin[x] y0,q[0]==0,q'[0]==0},q[x],x][[1]];
y1 = q[x]/.s1

s2 = DSolve[{r''[x]==Sin[x] y1,r[0]==0,r'[0]==0},r[x],x][[1]];
y2 = r[x]/.s2

s3 = DSolve[{t''[x]==Sin[x] y2,t[0]==0,t'[0]==0},t[x],x][[1]];
y3 = t[x]/.s3

ypert = y0 + y1 + y2 + y3

Remove["Global`*"]
s=NDSolve[{y''[x]== Sin[x] y[x],y[0]==1,y'[0]==1},y[x],{x,0,5}]

y0 = 1+x;
y1 = 2+x-2 Cos[x]-Sin[x]-x Sin[x];
y2 = (1/24) (51+33 x-6 (x^2)-2 (x^3)-48 Cos[x]-3 Cos[2 x]-3 x Cos[2 x]-48 Sin[x]-24 x Sin[x]+9 Sin[2 x]);
y3 =(1/432) (2284+882 x-216 (x^2)-36 (x^3)-2187 Cos[x]+432 x Cos[x]+216 (x^2) Cos[x]-108 Cos[2 x]-54 x Cos[2 x]+11 Cos[3 x]-1593 Sin[x]-1269 x Sin[x]+108 (x^2) Sin[x]+36 (x^3) Sin[x]+162 Sin[2 x]+3 Sin[3 x]+3 x Sin[3 x]);
ypert = y0 + y1 + y2 + y3;

Plot[{Evaluate[y[x]/.s],y0+y1,y0+y1+y2,y0+y1+y2+y3},{x,0,10},PlotStyle->Automatic]