Perturbation Solution of
(x^2)+ϵ x -1  =0

Find  an approximate solution of
(x^2)+ϵ x -1 =0
as
x(ϵ)=(∑_(n=0)^∞) (Subscript[x, n])(ϵ^n)
The equation obviously has an analytical solution, so the purpose
here is to demonstrate the perturbation method on an equation with 
a known solution, so that we can see that the method works.

Off[General::spell1]

First, define the equation:

(* Needs["Notation`"]; 
Symbolize[(NotationTemplateTag[(Subscript[\[SelectionPlaceholder], \[Placeholder]])])] *)

qe=x^2+ϵ x-1==0;

Find the exact solutions, which comes out as a list of 2 replacment rules:

qex=Solve[qe,x]

Show how to access the element of a list (in case you do not know):

qex[[2]]

qex[[2,1]]

Check that these solutions really do satisfy qe:

qe/.qex[[1,1]]//Simplify

qe/.qex[[2,1]]//Simplify

Define functions x(ϵ)for the two exact, analytical solutions

xe2[ϵ_]=x/.qex[[1]]

xe1[ϵ_]=x/.qex[[2]]

Now prepare to find an approximate  perturbation solution.  Note
the O[ϵ] symbol is powerful in Mathematica, try changing {n,0,2} to {n,0,10) and see what happens...

xs=Sum[(Subscript[x, n])(ϵ^n),{n,0,2}]+O[ϵ]^3

Substitute the series solution for xinto the equation qe:

peq=qe/.x->xs

LogicalExpand is really cool, and works because of the presence of the O[ϵ] symbol :

eqns=LogicalExpand[peq]

The "logical and" might just as well be a list of equations that need to be true.  In fact the equation can be accessed as list, for example:

eqns[[1]]

roots=Solve[eqns,{(Subscript[x, 0]),(Subscript[x, 1]),(Subscript[x, 2])}]

So here are the two approximate solutions to qe:

xa1[ϵ_]=Normal[xs]/.roots[[2]]

xa2[ϵ_]=Normal[xs]/.roots[[1]]

Compare the exact (green) and approximate (red) solutions:

Plot[{xa1[ϵ],xe1[ϵ]},{ϵ,-3.,3.},PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0]}]

Plot[{xa2[ϵ],xe2[ϵ]},{ϵ,-3.,3.},PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0]}]