# Solve for y1(x), y2(x) by the method of Frobenius, 
# using maple assist for substitution and solving the recurrences.

# Problem. Solve x^2 y'' - x(2+x) y' + 2x y = 0

# Answers: y1 = x^3(1+x/4+x^2/20+x^3/120+x^4/840+x^5/6720 + ...)
#          y2 = x^0(12 +12x + 6x^2 + x^4/2 + x^5/10 + ...)

de:=x^2*diff(y(x),x,x)-x*(2+x)*diff(y(x),x)+2*x*y(x)=0;
dsolve({de},y(x),series);

z1:=sum(a(n)*x^(n+3),n=0..infinity); # Frobenius trial solution, a(0)=1
subs(y(x)=z1,de);
expand(%);
simplify(lhs(%));

Q1:=n->a(n)*(n^2+5*n+6-2*n-6); # Coeff of x^(n+3)
Q2:=n->a(n)*(-n-3+2); # Coeff of x^(n+4)
Q1(0); # Coeff of x^3 should be zero (its the indicial equation)
Q1(n)+Q2(n-1)=0; # recursion relation for sequence a[n], valid for n>=1

# Use maple rsolve() to solve the first recurrence relation
eq:=Q1(n)+Q2(n-1)=0;
ic:=a(0)=1;
aa:=unapply(rsolve({eq,ic},a),n);
seq(aa(i),i=0..10);

z2:=0*ln(x)*y1(x)+sum(b(n)*x^(n+0),n=0..infinity); # Change 0 to k if it fails
subs(y(x)=z2,de);
expand(%);

simplify(lhs(%));
R1:=n->(n*(n-1)-2*n)*b(n); # Coeff of x^n
R2:=n->(-n+2)*b(n); # Coeff of x^(n+1)

# Use maple rsolve() to solve the second recurrence relation
eq:=R1(n)+R2(n-1)=0;
ic:=b(0)=12; # To match the dsolve output. Normally choose b(0)=1
bb:=unapply(rsolve({eq,ic},b),n);
seq(bb(n),n=0..10);