A Regular Perturbation Expansion
for 
y'=1+(1+ϵ)(y^2)

Solution of y'(x)=1 +(1+ϵ)(y^2) with 
y(0)=1.
We seek a solution of form
y(x)=(Subscript[y, 0])(x)+(Subscript[ϵy, 1])(x)+(ϵ^2)(Subscript[y, 2])(x)
Logan 1.4 on page 50 (1st edition).

Here is the proposed series solution:

ys= Sum[(Subscript[y, n])[t](ϵ^n),{n,0,2}]+O[ϵ]^3

D[ys,t]

Here is the equation with the proposed perturbation solution:

peq=D[ys,t]==1+(1+ϵ)ys^2

With the O[ϵ] symbol, the LogicalExpand does a powerful step:

eqns=LogicalExpand[peq]

Here is another way to look at it:

eqa=TableForm[ReplaceAll[eqns,And->List]]

Here is the ODE for (Subscript[y, 0]):

eqns[[1]]

...and it's solution:

sol0=DSolve[{eqns[[1]],(Subscript[y, 0])[0]==1},(Subscript[y, 0])[t],t]

Here is the ODE for (Subscript[y, 1]), making use of the known solution for (Subscript[y, 0]):

eqn1=eqns[[2]]/.sol0[[1]]

Here is the solution for (Subscript[y, 1]):

sol1=DSolve[{eqn1,(Subscript[y, 1])[0]==0},(Subscript[y, 1])[t],t]

sol0[[1]]

eqns[[3]]

eqn2=eqns[[3]]/.sol0[[1]]/.sol1[[1]]

sol2=DSolve[{eqn2,(Subscript[y, 2])[0]==0},(Subscript[y, 2])[t],t]

Now convert these replacement rules into approximate solutions, of increasing O[ϵ],
that we can plot:

ym0=ys+O[ϵ]^1/.sol0[[1]]

ym1=ys+O[ϵ]^2/.sol0[[1]]/.sol1[[1]]

ym2=ys+O[ϵ]^3/.sol0[[1]]/.sol1[[1]]/.sol2[[1]]

ya0[t_]=Normal[ym0]

ya1[t_,ϵ_]=Normal[ym1]

ya2[t_,ϵ_]=Normal[ym2]

Test these functions:

{ya0[.6],ya1[.6,.1],ya2[.6,.1]}

Lucky for us, we can find the exact solution:

ye=DSolve[{y'[t]==1+(1+ϵ)y[t]^2,y[0]==1},y[t],t]

and here it is as a convenient function:

yexact[t_,ϵ_]=y[t]/.ye[[1]]

Let's test it:

ya0[.6,.1]

ep=.1;

Here we plot the exact solution (green), the (Subscript[y, 0])(x)(red), (Subscript[y, 0])(x)+ϵ(Subscript[y, 1])(x)(blue),  and (Subscript[y, 0])(x) +(Subscript[ϵy, 1])(x)+(ϵ^2)(Subscript[y, 2])(x) (black).

Plot[{yexact[t,ep],ya0[t],ya1[t,ep],ya2[t,ep]},{t,0,.7},PlotStyle->{RGBColor[0,1,0],
RGBColor[1,0,0],RGBColor[0,0,1],RGBColor[0,0,0]}]