pfun=ParametricNDSolveValue[{y''[t]+y[t]==a Sin[2t ]+ Cos[2y[t]],y[0]==1,y'[0]==0},y,{t,0,10},{a}];
Show[Plot[Evaluate@Table[pfun[a][t],{a,-.6,.6,.1}],{t,0,5},PlotStyle->Lighter[Gray,0.3],PlotRange->All],Plot[pfun[0][t],{t,0,5},PlotStyle->{Lighter[Black,.4],Thick}],ImageSize->Medium]

Plot[Evaluate[pfun[0][t]+{-.5,0,.5}pfun'[0][t]],{t,0,5},PlotStyle->{{Dashed,Lighter[Blue,.4]},{Thick,Lighter[Black,.4]},{Dashed,Lighter[Red,.4]}},Filling->{1->{2},2->{3}},ImageSize->Medium]