Series Solutions

Mathematica commands used:
Sum, O, D, LogicalExpand, Solve, Collect, Table
Coefficient, Denominator, Factor Integer, Factorial

Series solutions near an ordinary point for a second order, linear, homogenous DE
    y'' + y = 0

3 terms in each of 2 linearly independent solns

First state the form of the solution to a desired number of terms.

Clear[y,a,y1,y2,i,x,lhs, ysoln]

y=Sum[a[i] x^i, {i,0,6}] + O[x]^7

Calculate the necessary derivatives

y1=D[y,x]

y2=D[y1,x]

Substitute y and its derivatives into the differential equation.

lhs=y2+y

Equate the coefficients for corresponding powers of x.

LogicalExpand[lhs==0]

Since this is second order, we should be able to solve for all but the first two terms in the series solution.

Solve[%,{a[2],a[3],a[4],a[5],a[6]}]

Substitute these back into y.

{ysoln}=y/.%

Collect the terms for the two linearly independent series solutions.

ysoln=Collect[ysoln,{a[0],a[1]}]

10 terms in each solution

Of course, this same technique is used if we wish to take these to 10 terms each.

Clear[y,a,y1,y2,i,x,lhs, ysoln, eqnset, varset, u1, u2]

y=Sum[a[i] x^i, {i,0,20}] + O[x]^21

y1=D[y,x]

y2=D[y1,x]

lhs=y2+y

eqnset=LogicalExpand[lhs==0]

varset=Table[a[i],{i,2,20}]

Solve[eqnset,varset]

{ysoln}=y/.%

ysoln=Collect[ysoln,{a[0],a[1]}]

u1=Coefficient[ysoln,a[0]]

u2=Coefficient[ysoln,a[1]]

Looking for patterns in the coefficients

Let's look at u1.  An obvious pattern is (-1) * something x^(2n).  Maybe factoring the coefficients will show something.

t=Table[Coefficient[u1,x^(2n)],{n,0,10}]

dlist=Denominator[t]

One thing you can try is getting prime factors.  In the following, we see that 24 = 2^3 * 3; 720=2^4 * 3^2 * 5, etc.

flist=FactorInteger[dlist]

You might try successive quotients.

flist2=Table[dlist[[i+1]]/dlist[[i]],{i,1,10}]

FactorInteger[%]

And you can always look for factorials.

Table[Factorial[n],{n,2,20}]

From this one can see that the even powers of x get even factorials and the odd powers get odd factorials.  This can be double checked.

v1=Sum[(-1)^n/Factorial[2n] x^(2n),{n,0,10}]

v2=Sum[(-1)^n/Factorial[2n+1] x^(2n+1),{n,0,10}]

u1-v1

u2-v2