(* ::Package:: *)

NSolve[x^5 - 16*x + 1 == 0, x]


eqn = x^5 - 16*x + \[Epsilon]; (* lhs of the equation to solve *)


ord = 3;              (* order of the small parameter to keep *)


X = Sum[a[n]*\[Epsilon]^n, {n,0,ord}];  (* expression for the root *)


ser = Series[eqn /. x->X, {\[Epsilon],0,ord}];  (* substitute the expression and expand *)


li = CoefficientList[ser, \[Epsilon]]  (* extrac the coefficients of the expansion *)


sol = Solve[li==0] (* equate all coefficients to 0 and solve the resulting equations *)


r = X /. sol (* substitute back *)


roots = r /. \[Epsilon]->1  (* recover the root of the original equation *)


N[roots] (*numerical values *)
