> restart;
# Example in Section 5.4
> xdot:= 2*x+y^2;
> ydot:= -2*y + x^2+y^2;
# 
# Picard Iteration Method
# Construction of the operator T
> Tx:= (x,y)-> unapply(-exp(2*t) *int(exp(-2*s)*y^2,s=t..infinity,0),t);
> Ty:= (x,y) -> unapply(exp(-2*t)*(sigma + int(exp(2*s)*(x^2+y^2),s=0..t)),t);
# Begin with the trivial function (0,0) and iterate
>  x1:=Tx(0,0);
> y1:=Ty(0,0);
> x2:=Tx(x1(s),y1(s)); 
> y2:=Ty(x1(s),y1(s));
> x3:= Tx(x2(s),y2(s));
> y3:= Ty(x2(s),y2(s));
# The graph of the stable manifold can be obtained by plotting our approximate
# solution for any value of t. Try t = 0--it is the simplest.
> x3(0);
> y3(0);
> plot([x3(0),y3(0), sigma = -1..1]);
> x4:= Tx(x3(s),y3(s)):
> y4:= Ty(x3(s),y3(s)):
> x4(0);
> y4(0);
> plot([x4(0),y4(0),sigma=-1..1]);
# 
# Power Series Method
# Now do this with power series method
# We assume that the stable manifold is given as x = h(y) = power series tangent to y axis
> h:= a*y^2+b*y^3+c*y^4+d*y^5;
# Invariance equation xdot(h(y),y) =  Dh(y) ydot(h(y),y)
> eq:= subs(x=h,xdot) -diff(h,y)* subs(x=h,ydot);
> e1:=coeff(eq,y^2);

> e2:=coeff(eq,y^3);

> e3:=coeff(eq,y^4);
> e4:=coeff(eq,y^5);
> solve({e1,e2,e3,e4},{a,b,c,d});
> subs(%,h);
> 
# 
#