# Residual computation for a system of two equations
restart; #top down execution should have a clean state to begin.
macro(e = epsilon); #saves typing "epsilon" every time.
Order := 4;
z := Vector(2,[3/5,4/5]); # z_0 = u_0
F := u -> Vector(2,
[ u[1]^2 + u[2]^2 - 1 - e*u[1]*u[2],
25*u[1]*u[2] - 12 + 2*e*u[1] ]);
Zer := F( [ x(e), y(e) ] ); #This asks for F evaluated at functions x(e)
# and y(e) that are yet unspecified.
diffeqs := { diff( Zer[1], e ), diff( Zer[2], e ) }; #This creates a set
# of two differential equations, one from each component of F.
# Each equation will contain both dx/de and dy/de.
iniconds := { x(0) = z[1] , y(0) = z(2) };
sol := dsolve( diffeqs union iniconds, {x(e), y(e)}, type=series );
Delta := eval( F( [x(e), y(e)] ), map(convert, sol, polynom) ):
map( series, Delta, e, Order+2 );