Chapter 5.

Series Solutions of ODEs 1

§5.1 Power Series Method
Power series (in powers of x − x0)
∞
∑
am(x − x0)m = a0 + a1(x − x0) + a2(x − x0)2 + · · · .
m=0

Examples(Maclaurin series)
1
= 1 + x + x2 + x3 + · · · (|x| < 1, geometric series)
1−x
x2 x3
e x
= 1+x+
+ + · · · (|x| < ∞)
2! 3!
x2 x4 x6
cos x = 1 − + − + · · · (|x| < ∞)
2! 4! 6!
x3 x5 x7
sin x = x − + − + · · · (|x| < ∞)
3! 5! 7!
Note
∞
∑
y= amxm = a0 + a1x + a2x2 + · · · (|x| < R)
m=0
∞
∑
y′ = mamxm−1 = a1 + 2a2x + 3a3x2 + · · · (|x| < R)
m=1
∞
∑
y ′′ = m(m − 1)amxm−2 = 2a2 + 3 · 2a3x + · · · (|x| < R)
m=2

Power series method

For a given ODE
y ′′ + p(x)y ′ + q(x)y = 0,
if we can represent p(x) and q(x) by power series in powers of x or x−x0,
then try
∞
∑ ∞
∑
y(x) = amx m
or y(x) = am(x − x0)m.
m=0 m=0

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 2

Example 2.
Solve the ODE by power series
y ′ − y = 0.
Sol. Set
∑∞
y= amxm = a0 + a1x + a2x2 + a3x3 + · · · .
m=0

Inserting the power series into the ODE, we obtain

a1 + 2a2x + 3a3x2 + 4a4x3 + 5a5x4 + · · ·
−(a0 + a1x + a2x2 + a3x3 + a4x4 + a5x5 + · · · ) = 0
and
a1 − a0 = 0, 2a2 − a1 = 0, 3a3 − a2 = 0, · · · .
We may express each coefficients a1, a2, · · · in terms of a0, which re-
mains arbitrary:
a1 a0 a2 a0
a1 = a0 , a 2 = = , a3 = = ,··· .
2 2! 3 3!
The solution is
x2 x3
y = a0(1 + x + + + · · · ) = a0ex.
2! 3!
More rapidly, inserting the power series into the ODE, we obtain

∞
∑ ∞
∑ ∞
∑
mamxm−1 − amxm = 0 or (mam − am−1)xm−1 = 0.
m=1 m=0 m=1

The recursion formula mam = am−1, m = 1, 2, · · · .

Thus, we have am = m! a0
, m = 0, 1, 2, · · · , and the solution is
∞
∑ ∞
∑ a0
y= amxm = xm = a0ex.
m=0 m=0
m!

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 3

Example 3. (A special Legendre equation)

Solve the ODE by power series
(1 − x2)y ′′ − 2xy + 2y = 0.

Sol. Set ∞
∑
y= amxm.
m=0
Inserting the power series into the ODE, we obtain
y ′′ = 2a2+ 6a3x+ 12a4x2+ 20a5x3+ 30a6x4 + · · ·
−x2y ′′ = − 2a2x2− 6a3x3− 12a4x4 − · · ·
−2xy ′ = − 2a1x− 4a2x2− 6a3x3− 8a4x4 − · · ·
2y = 2a0+ 2a1x+ 2a2x2+ 2a3x3+ 2a4x4 + · · ·
Adding terms of like powers of x, we have
x0 : a2 = −a0
x1 : a3 = 0
x2 : 12a4 = 4a2, a4 = 12 4
a2 = − 31 a0
x3 : a5 = 0 ( 1)
x4 : 30a6 = 18a4, a6 = 30 a4 = 30 − 3 a0 = − 15 a0
18 18

This gives the solution

x4 x6
y = a1x + a0(1 − x − 2
− − · · · ).
3 5

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 4

[Theory of the Power Series Method]

power series in powers of x − x0.
∑∞
am(x − x0)m = a0 + a1(x − x0) + a2(x − x0)2 + · · · (1)
m=0
∑
n ∞
∑
= am(x − x0) + m
am(x − x0)m
|m=0 {z } |m=n+1 {z }
sn (x) Rn (x)
where x0: center, am: coefficients, sn: n-th partial sum, Rn: remainder.
If lim sn(x1) = s(x1), then (1) is called convergent at x = x1.
n→∞
Otherwise, (1) is called divergent at x = x1.

[Convergence interval. Radius of convergence]

There are three cases:
Case 1: (1) may converges only at x = x0. [(1) always converges at x0.]
Case 2: (1) converges for all x such that |x − x0| < R
and diverges for all x such that |x − x0| > R for some R > 0.
In this case, R: radius of convergence of (1),
set of the values of x for which (1) converges: convergence interval.
am+1 −1
R = ( lim |am|1/m)−1 or R = ( lim | |)
m→∞ m→∞ am
if the limits exist.
Case 3: (1) converges for all x.

Fig. 105. Convergence interval of a power series with center x0

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 5

Example 4.
∑∞ ∞
∑ ∑∞
x xm 1 m
(a) e = (b) = x (c) m!xm
m=0
m! 1−x m=0 m=0
1/(m + 1)! 1
Sol. (a) lim | | = lim = 0, or R = ∞.
m→∞ 1/m! m→∞ m + 1
The power series converges for all x.

1
(b) lim | | = 1, or R = 1.
m→∞ 1
The power series converges only at |x| < 1.

(m + 1)!
(c) lim | | = ∞, or R = 0.
m→∞ m!
The power series converges only at x = 0.

[Existence of Power Series Solutions of ODEs]

Definition (Real Analytic Function)
A real function f (x) is analytic at x = x0 if it can be represented by a
power series in power of x − x0 with radius of convergence R > 0.

Theorem 1. (Existence of Power Series Solution)

y ′′ + p(x)y ′ + q(x)y = r(x)

If p, q, & r are analytic at x = x0, then every solution of the ODE is
analytic at x = x0 and can be represented by a power series in power of
x − x0 with radius of convergence R > 0.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 6

[Operations on Power Series]

Termwise Differentiation
∞
∑
y(x) = am(x − x0)m, |x − x0| < R
m=0
∞
∑
y(x)′ = mam(x − x0)m−1, |x − x0| < R
m=1
∞
∑
y(x)′′ = m(m − 1)am(x − x0)m−2, |x − x0| < R.
m=2

Termwise Addition and Multiplication

∞
∑
f (x) = am(x − x0)m, |x − x0| < R1,
m=0
∞
∑
g(x) = bm(x − x0)m, |x − x0| < R2.
m=0

∞
∑
=⇒ f (x) + g(x) = (am + bm)(x − x0)m,
m=0
|x − x0| < R = min{R1, R2}
∑∞
f (x)g(x) = (a0bm + a1bm−1 + · · · + amb0)(x − x0)m,
m=0
|x − x0| < R = min{R1, R2}.
Vanishing of all coefficients
∞
∑
am(x − x0)m = 0, |x − x0| < R
m=0

=⇒ a0 = a1 = · · · = 0.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 7

§5.2 Legendre’s Equation. Legendre Polynomials Pn(x)

[Legendre’s Equation]
(1 − x2)y ′′ − 2xy ′ + n(n + 1)y = 0.
Legendre’s equation can be rewritten as
′′ 2x ′ n(n + 1)
y − y + y = 0.
1 − x2 1 − x2
n(n+1)
Since − 1−x
2x
2 and 1−x2 are analytic at x = 0, Legendre’s equation has
power series solutions of the form
∞
∑
y= amxm.
m=0

Substituting y, y ′ and y ′′ into the ODE, we obtain

∞
∑ ∞
∑
(1 − x2) m(m − 1)amxm−2 − 2x mamxm−1
m=2 m=1
∞
∑
+ n(n + 1) amxm = 0
m=0

and
∞
∑ ∞
∑
m(m − 1)amx m−2
− m(m − 1)amxm
m=2 m=2
∞
∑ ∞
∑
− 2mamxm + n(n + 1) amxm = 0.
m=1 m=0

Setting m − 2 = s on the first series and m = s on the other series,

∞
∑ ∞
∑
(s + 2)(s + 1)as+2x − s
s(s − 1)asxs
s=0 s=2
∞
∑ ∞
∑
− 2sasxs + n(n + 1)asxs = 0.
s=1 s=0

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 8

The sum of the coefficients of each power of x on the left must be zero.
Thus
2a2 + n(n + 1)a0 = 0,
3 · 2a3 + [−2 + n(n + 1)]a1 = 0,
(s + 2)(s + 1)as+2 − [s(s − 1) + 2s − n(n + 1)]as = 0.
recurrence relation
(n − s)(n + s + 1)
as+2 = − as (s = 0, 1, 2, · · · )
(s + 2)(s + 1)

a2 = − n(n+1)
2!
a 0 , a 4 = (n−2)n(n+1)(n+3)
4!
a0, · · ·
a3 = − (n−1)(n+2)
3!
a 1 , a 5 = (n−3)(n−1)(n+2)(n+4)
5!
a1, · · ·

∴ y = a0y1 + a1y2
where
n(n + 1) (n − 2)n(n + 1)(n + 3)
y1 = 1 − x2 + x4 − + · · ·
2! 4!
(n − 1)(n + 2)
y2 = x − x3
3!
(n − 3)(n − 1)(n + 2)(n + 4)
+ x5 − + · · · .
5!
Note 1. These series converges for |x| < 1.
2. Since y1 and y2 are not proportional, they are linearly independent solu-
tions. y = a0y1 + a1y2 is a general solution of the Legendre’s equation
on −1 < x < 1.
3. n: even, y1(x) is a polynomial of degree n,
n: odd, y2(x) is a polynomial of degree n.
These polynomials, multiplied by some constants, are called Legendre poly-
nomials.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 9

Let n be a positive integer.

Choose the coefficient an of the highest power xn as
(2n)! 1 · 3 · 5 · · · (2n − 1)
an = = .
2n(n!)2 n!
The choice makes Pn(1) = 1.
In the recurrence relation, putting s = n − 2
(s + 2)(s + 1)
as = − as+2, (s ≤ n − 2)
(n − s)(n + s + 1)
gives
n(n − 1)
an−2 = − an
2(2n − 1)
n(n − 1) (2n)!
= −
2(2n − 1) 2n(n!)2
n(n − 1) 2n(2n − 1)(2n − 2)!
= −
2(2n − 1) 2nn(n − 1)!n(n − 1)(n − 2)!
(2n − 2)!
= − n .
2 (n − 1)!(n − 2)!
Similarly, putting s = n − 4, s = n − 6, · · · , s = n − 2m, we have
(2n − 2m)!
an−2m = (−1)m .
2nm!(n − m)!(n − 2m)!
Legendre polynomials Pn(x) of degree n
∑
M
(2n − 2m)!
Pn(x) = (−1)m xn−2m
m=0
2nm!(n − m)!(n − 2m)!
(2n)! (2n − 2)!
= xn
− x n−2
+ ···
2n(n!)2 2n1!(n − 1)!(n − 2)!
where M = n/2 or (n − 1)/2.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 10

Legendre polynomials
P0(x) = 1, P1(x) = x
1 1
P2(x) = (3x − 1),
2
P3(x) = (5x3 − 3x)
2 2
1 1
P4(x) = (35x4 − 30x2 + 3), P5(x) = (63x5 − 70x3 + 15x)
8 8

··· , ···

Note Pn(x), n = 0, 1, 2, · · · are orthogonal, that is,

∫ 1
Pn(x)Pm(x)dx = 0, (n ̸= m).
−1

Fig. 107. Legendre polynomials

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 11

§5.3 Extended Power Series Method: Frobenius Method

Theorem 1.(Frobenius method)
Let b(x) and c(x) be analytic at x = 0. Then
b(x) c(x)
y ′′ + y′ + y=0
x x2
has at least one solution y represented in the form
∞
∑
y(x) = x r
amxm (a0 ̸= 0)
m=0
= x (a0 + a1x + a2x2 + · · · )
r

where r may be any (real or complex) number.

The ODE also has a second solution (linearly independent)
∑
(similar to xr ∞
m=0 am x
m
or may contain a logarithmic term).
Ex. Bessel’s equation
′′ 1 ′ x2 − ν 2
y +
y + y=0
x x2
Since b(x) = 1, c(x) = x2 − ν 2 are analytic at x0 = 0, Frobenius
method can be applied.

Definition (Regular and Singular Points)

A regular point of
y ′′ + p(x)y ′ + q(x)y = 0
is a point x0 at which p and q are analytic. Then the power series method
can be applied. If x0 is not regular, it is called singular.
Similarly, a regular point of the ODE
h̃(x)y ′′ + p̃(x)y ′ + q̃(x)y = 0
is a point x0 at which h̃, p̃ and q̃ are analytic and h̃(x0) ̸= 0.
If x0 is not regular, it is called singular.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 12

[Indicial Equation, Indicating the Form of Solutions]

Assume that b(x) and c(x) are analytic at x = 0.
b(x) c(x)
y ′′ + y′ + y=0
x x2

⇒ x2y ′′ + xb(x)y ′ + c(x)y = 0.

Since b(x) and c(x) are analytic at x = 0, we expand b(x) and c(x) in
power series,
b(x) = b0 + b1x + b2x2 + · · · ,
c(x) = c0 + c1x + c2x2 + · · · .
By the Frobenius method, the ODE has a solution represented in the form
∞
∑
y= amxm+r , a0 ̸= 0.
m=0

Differentiate y, then
∞
∑
y′ = (m + r)amxm+r−1
m=0
= xr−1
[ra0 + (r + 1)a1x + · · · ] ,
∑∞
y ′′ = (m + r)(m + r − 1)amxm+r−2
m=0
= xr−2
[r(r − 1)a0 + (r + 1)ra1x + · · · ] ,
by inserting all these power series into the ODE, we obtain
xr [r(r − 1)a0 + · · · ] + (b0 + b1x + · · · )xr (ra0 + · · · )
+ (c0 + c1x + · · · )xr (a0 + a1x + · · · )
= [r(r − 1)a0 + b0ra0 + c0a0]xr + · · · = 0.
Thus
[r(r − 1) + b0r + c0]a0 = 0.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 13

Since a0 ̸= 0, we have
r(r − 1) + b0r + c0 = 0 : indicial equation of the ODE
roots: r = r1, r2

Note Here, b0 = b(0), c0 = c(0).

Case 1 : Distinct roots, r1 − r2 is not an integer.

Case 2 : Double root. r1 = r2.
Case 3 : Distinct roots, r1 − r2 is an integer.

Theorem 2. (Frobenius Method. Basis of Solutions. Three Cases)

Suppose that the ODE
b(x) c(x)
y ′′ + y′ +
y=0
x x2
satisfies the assumptions in Theorem 1. Let r1 and r2 be the roots of the
indicial equation. Then
Case 1 (r1 − r2 is not an integer.): A basis is
y1(x) = xr1 (a0 + a1x + a2x2 + · · · ),
y2(x) = xr2 (A0 + A1x + A2x2 + · · · ).
Case 2 (r1 = r2 = r = 12 (1 − b0)) : A basis is
y1(x) = xr (a0 + a1x + a2x2 + · · · ),
y2(x) = y1(x) ln x + xr (A1x + A2x2 + · · · ),
where x > 0.
Case 3 (r1 − r2(> 0) is an integer) : A basis is
y1(x) = xr1 (a0 + a1x + a2x2 + · · · )
y2(x) = ky1(x) ln x + xr2 (A0 + A1x + A2x2 + · · · )
where x > 0 and k may turn out to be zero.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 14

[Recall : Method of Reduction of Order]

When y1 is a solution of the ODE
y ′′ + p(x)y ′ + q(x)y = 0,
a second linearly independent solution is
∫
1 −∫ p dx
y2 = y1 u = y1 2
e dx.
y1
[Appendix]
1−b0
Consider Case 2. That is, r = .
The indicial equation is
2
1
r(r − 1) + b0r + c0 = r 2 + (b0 − 1)r + (b0 − 1)2 = 0.
4
Let
y1(x) = xr (a0 + a1x + a2x2 + · · · ), a0 ̸= 0.
Set y2 = uy1 by the method of reduction of order and a0 = 1. Then we
have ∫
1 − ∫ b(x) dx
u = 2
e x dx
y1
∫
1 − ∫ b0+b1x+b2x2+··· dx
= 2
e x dx
y1
∫
1 − ∫ b0 +b1+b2x+···dx
= 2
e x dx
y1
∫
x−b0 e−(b1x+b2/2x +··· )
2

= dx
x (1 + a1x + a2x + · · · )
2r 2 2
∫ −(b0+2r)
x (1 + d1x + d2x2 + · · · )
= 2 + · · · )2
dx
∫ (1 + a 1 x + a 2 x
= x−1(1 + e1x + e2x2 + · · · )dx
e2
= ln x + A1x + A2x + · · · , 2
(A1 = e1, A2 = , · · · ).
2
Hence
y2 = uy1 = y1 ln x + xr (A1x + A2x2 + · · · ).

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 15

Example 1. Euler-Cauchy Equation

Solve the Euler-Cauchy Equation
x2y ′′ + b0xy ′ + c0y = 0.

Sol. y = xr ⇒ r(r − 1) + b0r + c0 = 0 : indicial equation

(i) r1 ̸= r2 ⇒ y1 = xr1 , y2 = xr2 : y = c1xr1 + c2xr2 .
(real or complex numbers)
(ii) r1 = r2 = r ⇒ y1 = xr , y2 = xr ln x : y = (c1 + c2 ln x)xr .
Example 2.
Solve the ODE
x(x − 1)y ′′ + (3x − 1)y ′ + y = 0.

Sol. Note that

3x−1 x
y ′′ + x−1
y′ + x−1
y = 0.
x x2
3x − 1 x
Since b(x) = and c(x) = are analytic at x = 0,
x−1 x−1
by the Frobenius method, set
∑∞
y= amxm+r , (a0 ̸= 0).
m=0
indicial equation
r(r − 1) + b0r + c0 = r(r − 1) − r = r 2 = 0, r = 0 (double root).
∞
∑
First Solution: Inserting y = amxm into the ODE, we obtain
m=0
∞
∑ ∞
∑
m(m − 1)amx m
− m(m − 1)amxm−1
m=2 m=2
∞
∑ ∞
∑ ∞
∑
+3 mamxm − mamxm−1 + amxm = 0.
m=1 m=1 m=0

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 16

Setting m = s + 1 on the second and fourth series and m = s on the

others,
∞
∑ ∞
∑ ∞
∑
s(s − 1)asxs − (s + 1)sas+1xs + 3 sasxs
s=2 s=1 s=1
∞
∑ ∞
∑
− (s + 1)as+1x + s
asxs
s=0 s=0
∞
∑
= (a0 − a1) + 4(a1 − a2)x + (s + 1)2[as − as+1]xs
s=2
= 0,
and so
a0 = a1 = a2 = a3 = · · · .
By choosing a0 = 1, we obtain
∞
∑ 1
y1(x) = xm = (|x| < 1).
m=0
1−x
u
Second Solution : Set y2 = uy1 = 1−x by the method of reduction of
order, we have
∫
1 − ∫ x(x−1)
3x−1
dx
u = 2
e dx
∫ 1y
∫ 1 2
2 − x
= (1 − x) e + x−1 dx
dx
∫ ∫
1 1
= (1 − x) [
2
]dx = dx
x(x − 1)2 x
= ln x (why!)
Hence
ln x
y2 = .
1−x
y1 and y2 are linearly independent and form a basis on 0 < x < 1 (as well
as on 1 < x < ∞). Thus the general solution is y = 1−x 1
(c1 + c2 ln x).

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 17

Example 3.
Solve the ODE
(x2 − x)y ′′ − xy ′ + y = 0.

Sol. Note that x x

y ′′ − x−1
y′ + x−1
y = 0.
x x2
x x
Since b(x) = − and c(x) = are analytic at x = 0,
x−1 x−1
by the Frobenius method, put
∑∞
y= amxm+r , (a0 ̸= 0).
m=0
indicial equation
r(r − 1) + b0r + c0 = r(r − 1) = 0, r1 = 1, r2 = 0.
∞
∑
First Solution : Substitution y = amxm+1, (r1 = 1) gives
m=0
∞
∑ ∞
∑
(m + 1)mamxm+1 − (m + 1)mamxm
m=1 m=1
∞
∑ ∑∞
− (m + 1)amxm+1 + amxm+1 = 0.
m=0 m=0
Setting m = s + 1 on the second series and m = s on the others,
∞
∑ ∑∞
(s + 1)sasx s+1
− (s + 2)(s + 1)as+1xs+1
s=1 s=0
∞
∑ ∞
∑
− (s + 1)asxs+1 + asxs+1
s=0 s=0
∞
∑
= −2a1x + [s2as − (s + 2)(s + 1)as+1]xs+1
s=1
= 0.

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 18

Thus a1 = a2 = a3 = · · · = 0 and choosing a0 = 1, we have

y1(x) = x.
Second Solution : Set y2 = uy1 = ux by the method of reduction of
order, we have
∫ ∫ −1
1 −
u = e x−1 dx dx
y12
∫
x−1
= dx
x2
1
= ln x + .
x

∴ y2 = ux = x ln x + 1.
Since y1 and y2 are linearly independent, they form a basis. The
general solution is
y = c1x + c2(x ln x + 1).

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 19

§5.4 Bessel’s Equation. Bessel Functions Jν (x)

[Bessel’s Equation]
x2y ′′ + xy ′ + (x2 − ν 2)y = 0 (ν ≥ 0)

′′ 1 ′ x2 − ν 2
=⇒ y + y + y = 0.
x x2
Since b(x) = 1 and c(x) = x2 − ν 2 are analytic at x = 0, by the
Frobenius method, the ODE[Bessel’s Equation] has a solution of the form
∞
∑
y= amxm+r , a0 ̸= 0.
m=0

Inserting the power series into the ODE, we obtain

∞
∑ ∞
∑
(m + r)(m + r − 1)amx m+r
+ (m + r)amxm+r
m=0 m=0
∞
∑ ∞
∑
+ amxm+r+2 − ν 2
amxm+r = 0
m=0 m=0

and (m = s in first, second, fourth series, m = s − 2 in the third series)

[r(r − 1)a0 + ra0 − ν 2a0]xr

+[(r + 1)ra1 + (r + 1)a1 − ν 2a1]xr+1

∞
∑ ∞
∑
+ (s + r)(s + r − 1)asxs+r + (s + r)amxs+r
s=2 s=2
∞
∑ ∞
∑
+ as−2x s+r
−ν 2
asxs+r = 0.
s=2 s=2

indicial equation r(r − 1) + r − ν 2 = r 2 − ν 2 = 0

r1 = ν, r2 = −ν (ν ≥ 0)
Computational Science & Engineering (CSE) T. Jeong
 Chapter 5. Series Solutions of ODEs 20

When r = r1 = ν,
(ν + 1)νa1 + (ν + 1)a1 − ν 2a1 = (2ν + 1)a1 = 0,
(s + 2ν)sas + as−2 = 0, s = 2, 3, · · · .
Since ν ≥ 0, a1 = 0 and
1
as = − as−2.
(s + 2ν)s
Thus a3 = a5 = · · · = 0, and for s = 2m,
1
a2m = − a2m−2, m = 1, 2, · · · .
2m(2m + 2ν)
This gives
a0
a2 = − ,
22(ν + 1)
a2 a0
a4 = − =− ,
22(ν + 2) 242!(ν + 1)(ν + 2)

··· ,

(−1)ma0
a2m = , m = 1, 2, · · · .
22mm!(ν + 1)(ν + 2) · · · (ν + m)

Hence we have a solution

∞
∑
Jν (x) ≡ y1(x) = a2mx2m+ν
m=0
∞
∑ (−1)m
= a0 x2m+ν .
m=0
22mm!(ν + 1)(ν + 2) · · · (ν + m)

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 21

[Bessel functions Jn(x) for integer ν = n]

1
Choosing a0 = 2n n!
(why!),
(−1)ma0
a2m =
22mm!(n + 1)(n + 2) · · · (n + m)
(−1)m
= 2m+n ,
2 m!(n + m)!
and
∑
Jn(x) = xn ∞ (−1)m x2m
m=0 22m+n m!(n+m)! .
By the ratio test, the series converges for all x. Jn(x) is called the Bessel
function of the first kind of order n.
Example 1. Bessel functions
Bessel functions of order 0 and 1 are given as follows:
∞
∑ (−1)mx2m
J0(x) =
m=0
22m(m!)2
x2 x4 x6
= 1− + 4 − 6 + −···
22(1!)2 2 (2!)2 2 (3!)2
∞
∑ (−1)mx2m+1
J1(x) =
m=0
22m+1m!(m + 1)!
x x3 x5 x7
= − + − + −··· .
2 231!2! 252!3! 273!4!

Fig. 110. Bessel functions of the first kind J0 and J1

√ ( )
Note Jn(x) ∼ 2
πx
cos x − nπ
2
− π
4

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 22

[Gamma Function]
∫ ∞
Γ(ν) = e−ttν−1 dt, ν > 0
0
Note By integrating by parts we obtain
Γ(ν + 1) = νΓ(ν).
Note that Γ(n + 1) = nΓ(n) = n! (n = 1, 2, · · · ).
1
Choosing a0 = 2ν Γ(ν+1) ,
(−1)ma0
a2m =
22mm!(ν + 1)(ν + 2) · · · (ν + m)
(−1)m
=
22mm!(ν + 1)(ν + 2) · · · (ν + m)2ν Γ(ν + 1)
(−1)m
= 2m+ν ,
2 m!Γ(ν + m + 1)
and ∞
∑ (−1)mx2m
ν
Jν (x) = x .
m=0
22m+ν m!Γ(ν + m + 1)
By the ratio test, the series converges for all x. Jν (x) is called the Bessel
function of the first kind of order ν.

[Properties from series]

Theorem 1. (Derivatives, Recursions)
Derivatives of Jν (x) with respect to x
[a] [xν Jν (x)]′ = xν Jν−1(x)
[b] [x−ν Jν (x)]′ = −x−ν Jν+1(x)
Recurrence relations
2ν
[c] Jν−1(x) + Jν+1(x) = Jν (x)
x
[d] Jν−1(x) − Jν+1(x) = 2Jν′ (x)

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 23

Proof. (a) Note that

∞
∑
ν (−1)mx2m+2ν
x Jν (x) = .
m=0
22m+ν m!Γ(ν + m + 1)
Differentiating and using Γ(ν + m + 1) = (ν + m)Γ(ν + m), we have
∞
∑ (−1)m2(m + ν)x2m+2ν−1
(xν Jν (x))′ =
m=0
22m+ν m!Γ(ν + m + 1)
∑∞
(−1)mx2m+2ν−1
=
m=0
22m+ν−1m!Γ(ν + m)
∞
∑ (−1)mx2m
= xν xν−1
m=0
22m+ν−1m!Γ(ν + m)
= xν Jν−1(x).
(b) Note that
∞
∑
−ν (−1)mx2m
x Jν (x) =
m=0
22m+ν m!Γ(ν + m + 1)
and the first term is constant. Differentiating and setting m = s + 1, we
have
∑∞
−ν ′ (−1)m2mx2m−1
(x Jν (x)) =
m=1
22m+ν m!Γ(ν + m + 1)
∑∞
(−1)mx2m−1
=
m=1
22m+ν−1(m − 1)!Γ(ν + m + 1)
∑∞
(−1)s+1x2s+1
=
s=0
22s+ν+1s!Γ(ν + s + 2)
∑∞
−ν ν+1 (−1)sx2s
= −x x
s=0
22s+ν+1s!Γ(ν + s + 2)
= −x−ν Jν+1(x).

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 24

(c),(d) We perform the differentiation in (a), and do the same in (b) and
multiply the result on both sides by x2ν . This gives
(a∗) νxν−1Jν + xν Jν′ = xν Jν−1,
(b∗) − νxν−1Jν + xν Jν′ = −xν Jν+1.
Subtracting (b∗) from (a∗) and dividing the result by xν gives (c).
Adding (a∗) and (b∗) and dividing the result by xν gives (d).

Example 2. Application
1. Consider the recurrence formula[Theorem 1 (c)]
2ν
Jν+1(x) = Jν (x) − Jν−1(x). (I)
x
For ν = 1, J2(x) = 2J1(x)/x − J0(x), so that J2 can be obtained
from tables of J0 and J1.

2. Theorem 1 (b) gives

∫ 2 2
−3 −3
I = x J4(x)dx = −x J3(x)
1 1
1
= − J3(2) + J3(1).
8
From (I) we get J3 = 4x−1J2 − J1 and J2 = 2x−1J1 − J0. Together,
J3 = (8x−2 − 1)J1 − 4x−1J0.
Hence,
1
I = − J3(2) + J3(1)
8
1
= − (J1(2) − 2J0(2)) + 7J1(1) − 4J0(1).
8

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 25

Example 3. Bessel function Jν with ν = ± 12 , ± 23 , · · ·

Bessel functions Jν of orders ± 12 , ± 32 , · · · can be expressed by finitely
many cosines and sines and powers of x. In particular,
√ √
2 2
[a] J1/2(x) = sin x, [b] J−1/2(x) = cos x.
πx πx

Proof. (a) Note that

∞ √ √ ∞
∑ (−1)mx2m x 2 ∑ (−1)mx2m+1
J1/2(x) = 2m+1/2 m!Γ(m + 3 )
= .
m=0
2 2
x m=0 22m+1m!Γ(m + 32 )
Using
2mm! = 2m(2m − 2)(2m − 4) · · · 4 · 2,
3 1 1 3 1 1
2m+1Γ(m + ) = 2m+1(m + )(m − ) · · · · Γ( )
2 2 2 2 √2 2
= (2m + 1)(2m − 1) · · · 3 · 1 · π,
√
(here we have used Γ( 12 ) = π.) we have
√ ∞ √
2 ∑ (−1)mx2m+1 2
J1/2(x) = = sin x.
πx m=0 (2m + 1)! πx
(b) By Theorem 3(a), we have
√
√ 2
x1/2J−1/2(x) = [ xJ1/2(x)]′ = cos x.
π
This proves (b). From Theorem 3 and (a), (b) follow further formulas suc-
cessively.
Note
√
1 2 sin x
J 3 (x) = J 1 (x) − J− 1 (x) = ( − cos x),
2 x 2 2 πx√ x
1 2 cos x
J− 3 (x) = − J− 1 (x) − J 1 (x) = − ( + sin x).
2 x 2 2 πx x

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 26

[General Solution. Linear Dependence]

Note 1. Bessel’s Equation : x2y ′′ + xy ′ + (x2 − ν 2)y = 0.
2. Indicial equation : r 2 − ν 2 = 0.
3. Bessel function Jν (x) of the first kind of order ν is one of the solutions.
4. For a general solution, we need a second linearly independent solution.
∑∞
−ν (−1)mx2m
J−ν (x) = x
m=0
22m−ν m!Γ(m − ν + 1)
5. If ν is not an integer, Jν (x) and J−ν (x) are linearly independent
because the first term in Jν (x) and the first term in J−ν (x) are finite
nonzero multiples of xν and x−ν , respectively. Thus, if ν is not an integer,
a general solution of Bessel’s equation for all x ̸= 0 is
y(x) = c1Jν (x) + c2J−ν (x).

6. If ν = n is an integer, then c1Jn(x) + c2J−n(x) is not a general

solution because of linear dependence.
7. Linear dependence : J−n(x) = (−1)nJn(x) n = 1, 2, · · · .
Proof. Notice that Γ(m) is infinite for m = −1, −2, · · · .
Let n be a positive integer. Since Γ(m − n + 1) is infinite for m =
0, 1, · · · , n − 1, we obtain by using Γ(m − n + 1) = (m − n)!
∑∞
(−1)mx2m−n
J−n(x) =
m=0
22m−nm!Γ(m − n + 1)
∑∞
(−1)mx2m−n
=
m=n
22m−nm!(m − n)!
∑ ∞
(−1)n+sx2s+n
= 2s+n
(m = n + s)
s=0
2 (n + s)!s!
∑∞
n n (−1)sx2s
= (−1) x
s=0
22s+n(n + s)!s!
= (−1)nJn(x).

Computational Science & Engineering (CSE) T. Jeong

 Chapter 5. Series Solutions of ODEs 27

§5.5 Bessel functions of the second kind Yν (x)

(Brief discussion)
[n = 0: Bessel function of the second kind Y0(x)]
When n = 0, Bessel equation is
xy ′′ + y ′ + xy = 0.
indicial equation r(r − 1) + r = r 2 = 0, r = 0(double root)
We have one solution
∑∞
(−1)mx2m
y1(x) = J0(x) = 2m 2
.
m=0
2 (m!)
From Theorem 2 in Sec. 5.3, the second solution should be of the form
∑∞
n
y2(x) = J0(x) ln x + x Amxm.
m=1
Substituting y2 into the ODE, we have
A1 = A3 = · · · = 0,
(−1)m−1hm
A2m = , m = 1, 2, · · · ,
22m(m!)2
where h1 = 1, hm = 1 + 1
2
+ ··· +1
m
, and
∞
∑ (−1)m−1hm
y2(x) = J0(x) ln x + x2m.
m=1
22m(m!)2
Since J0 and y2 are linearly independent, they form a basis of the ODE for
x > 0. If we replace y2 by an independent particular solution of the form
a(y2 + bJ0)
where a(̸= 0), b are constants, we obtain another basis of the ODE.
Choose
2
a= and b = γ − ln 2
π
where γ = 0.577 · · · is the limit of 1 + 12 + · · · + 1s − ln s, called the
Euler constant.
Computational Science & Engineering (CSE) T. Jeong
 Chapter 5. Series Solutions of ODEs 28

The particular solution is[

2
( x ) ∑∞ (−1)m−1hm 2m]
Y0(x) = π J0(x) ln 2 + γ + m=1 22m(m!)2 x ,
called the Bessel function of the second kind of order zero or Neu-
mann’s function of order zero.
For ν = n = 1, 2, · · · a second solution can be obtained by manipulations
similar to those for n = 0.
We introduce a standard second solution Yν (x) defined for all ν by
1 [ ]
Yν (x) = Jν (x) cos νπ − J−ν (x) ,
sin νπ

Yn(x) = lim Yν (x).

ν→n

This function is called the Bessel function of the second kind of orderν
or Neumann’s function of order ν.
Notebthat Tν and Yν are linearly independent for all ν(and x > 0), and
for n = 0, 1, 2, · · ·
∞
( x
2 ) xn ∑ (−1)m−1(hm + hm+n) 2m
Yn(x) = Jn(x) ln + γ + x
π 2 π m=0 22m+nm!(m + n)!
x−n ∑ (n − m − 1)!
n−1
− x2m.
π m=0
22m−nm!

Theorem 1. General solution of Bessel’s equation

A general solution of Bessel’s equation for all ν (and x > 0) is
y(x) = C1Jν (x) + C2Yν (x).

Example
Solve the following ODEs.
(a) x2y ′′ +xy ′ +(x2 − 19 )y = 0 (b) x2y ′′ +xy ′ +(x2 −4)y = 0
Sol. (a) y(x) = c1J 1 (x) + c2J− 1 (x) or y(x) = c1J 1 (x) + c2Y 1 (x)
3 3 3 3
(b) y(x) = c1J2(x) + c2Y2(x)

Computational Science & Engineering (CSE) T. Jeong