# Perturbation solution of F(u;epsilon) = 0
restart; #top down execution should have a clean state to begin.
macro(e = epsilon); #saves typing "epsilon" every time.
F := z -> z^5-e*z-1;
# Zeroth order solution, by inspection, is
z := 1; #solve(eval(F(z),e=0),z);
A := coeff(series((D(F))(z), e, 1), e, 0); #A must not be 0 for regularity
N := 3; #number of terms
Delta := F(z); #initial residual, must be O(e)
# Now, the iteration:
for k to N do
u := -coeff(series(Delta, e, k+1), e, k);
z := z + u*e^k/A;
Delta := F(z);
end do;
z;
series(Delta, e, N+3);