Finding terms of the perturbation solution

In Version 12 you can now do

AsymptoticDSolveValue[{y'[t] == 1 + (1 + e) y[t]^2, y[0] == 1}, y[t], t, {e, 0, 4}]

Note that the result for this problem uses complex numbers so when you plot the results you want to use Chop[] or Re[] when converting the exact solution to a numerical one. MMA returns answers with tiny (10^-16) or even 0.0 complex parts.

You may start with something like that

equ = {-y'[t] + 1 + (1 + \[Epsilon]) y[t]^2};
y[t_] := Sum[x[i][t] \[Epsilon]^i, {i, 0, 10}] // Evaluate;

First order solution (use SeriesCoefficient function to expand equation with series substituted)

x[0] = x[0] /. 
  First@DSolve[{First@SeriesCoefficient[equ, {\[Epsilon], 0, 0}] == 0,
      x[0][0] == 1}, x[0], t];

and then loop over higher order terms (all with zero BC)

Do[
  x[k] = x[k] /. 
    First@DSolve[{First@SeriesCoefficient[equ, {\[Epsilon], 0, k}] == 
        0, x[k][0] == 0}, x[k], t], {k, 1, 4}];

Set not defined coefficients to zero by

x[_] = Function[t, 0];

Finally, compare with explicit solution, either by plotting

Block[{\[Epsilon] = 0.1}, 
 Plot[{y[t], 
   Tan[t Sqrt[1 + \[Epsilon]] + ArcTan[Sqrt[1 + \[Epsilon]]]]/Sqrt[
   1 + \[Epsilon]]}, {t, -3, 1}, PlotRange -> Automatic,
 PlotLegends -> "Expressions"]]

or comparing series expansions

Equal @@ Normal@
   Series[{y[t], 
     Tan[t Sqrt[1 + \[Epsilon]] + ArcTan[Sqrt[1 + \[Epsilon]]]]/Sqrt[
     1 + \[Epsilon]]}, {\[Epsilon], 0, 4}] // Simplify

True