Analysis and object oriented implementation of the Kovacic algorithm

Analysis and object oriented implementation
of the Kovacic algorithm

September 9, 2023

Abstract

This paper gives a detailed overview and a number of worked out examples illustrating the Kovacic [1] algorithm for solving second order linear diﬀerential equation $A (x) y^{″} + B (x) y^{'} + C (x) y = 0$ where $A, B, C$ are rational functions with complex coeﬃcients in the independent variable $x$ . All three cases of the algorithm were implemented in a software package based on an object oriented design and complete source code listing given in the appendix with usage examples. Implementation used the Maple computer algebra language.¹. This package was then used to analyze the distribution of Kovacic algorithm cases on $3000$ diﬀerential equations.

1 Introduction
2 Description of algorithm for each case
  2.1 Case one
   2.1.1 step 1
   2.1.2 step 2
   2.1.3 step 3
   2.1.4 Algorithm ﬂow chart for case one
  2.2 Case two
   2.2.1 step 1
   2.2.2 step 2
   2.2.3 step 3
   2.2.4 Algorithm ﬂow chart for case two
  2.3 Case three
   2.3.1 step 1
   2.3.2 step 2
   2.3.3 step 3
   2.3.4 Algorithm ﬂow chart for case three
  2.4 Statistics and discussion of results obtained using Kovacic algorithm
3 Worked example for each case
  3.1 case one
   3.1.1 Example 1
   3.1.2 Example 2
   3.1.3 Example 3
  3.2 case two
   3.2.1 Example 1
   3.2.2 Example 2
  3.3 case three
   3.3.1 Example 1
   3.3.2 Example 2
4 Conclusion
5 Appendix
  5.1 Instructions and examples using the Kovacic package
  5.2 Source code

1 Introduction

Kovacic [1] gave an algorithm for ﬁnding a closed form Liouvillian² solution to any linear second order diﬀerential equation $A y^{″} + B y^{'} + C y = 0$ if such a solution exists. Smith [2] gave an implementation based on a modiﬁed version of Kovacic algorithm by Saunders [3].

The current implementation is based on the original paper by Kovacic and uses the new object oriented features in Maple. The accompanied software package have been tested on $3000$ diﬀerential equations with each solution veriﬁed using Maple’s odetest. The test suite is included as a separate module. The Appendix describes how to use the software.

The Kovacic algorithm ﬁnds one (basis) solution of $A y^{″} + B y^{'} + C y = 0$ . The second basis solution is found using reduction of order. The general solution is a linear combination of the two basis solutions found.

The algorithm starts by writing the input ode $A y^{″} + B y^{'} + C y = 0$ as $\begin{array}{r} (1) & y^{″} + a y^{'} + b y = 0 \end{array}$

Where $a = \frac{B}{A}, b = \frac{C}{A}$ . The substitution $\begin{aligned} (2) & z & = y e^{\frac{1}{2} \int a d x} \end{aligned}$

is then applied to (1) which transforms it to a second order ode in the new dependent variable $z (x)$ without the ﬁrst derivative $\begin{aligned} (3) & z^{″} & = r z \end{aligned}$

$r$ in the above is given by $\begin{aligned} (4) & r & = \frac{1}{4} a^{2} + \frac{1}{2} a^{'} - b \end{aligned}$

It is ode (3) which is solved by the algorithm and not (1). Equation (3) will be called the DE from now on.

If a solution $z (x)$ to the DE is found, then the ﬁrst basis solution to the original ode is obtained using the transformation (2) in reverse $\begin{aligned} y & = z e^{- \frac{1}{2} \int a d x} \end{aligned}$

The second solution is found using reduction of order.

These are the four possible cases to consider.

DE has solution $z = e^{\int ω d x}$ with $ω \in C (x)$ .
DE has solution $z = e^{\int ω d x}$ with $ω$ polynomial over $C (x)$ of degree $2$ .
Solutions of DE are algebraic over $C (x)$ .
DE has no Liouvillian solution.

Before describing how the algorithm works, there are necessary (but not suﬃcient) conditions that determine which case the DE satisﬁes. Only those cases that meet the necessary conditions will be attempted.

The following are the necessary conditions for each case. To check each case, let $r = \frac{s}{t}$ where $gcd (s, t) = 1$ . The order of $r$ at $\infty$ (from now on referred to as $O (\infty)$ ) is deﬁned as $\deg (t) - \deg (s)$ . The poles of $r$ and the order of each pole need to be determined.

Knowing the order of the poles of $r$ and $O (\infty)$ is all what is needed to determine the necessary conditions for each case. These conditions are the following

Case $1$ . Either no pole exists, or if a pole exists, the order must be either one or even. If $O (\infty)$ is less $3$ , then it must be even otherwise it can be even or odd.
Case $2$ . $r$ must have at least one pole either of order 2 or odd order greater than 2. There are no conditions on $O (\infty)$ .
Case $3$ . $r$ must have a pole either of order $1$ or $2$ . No other order is allowed. $O (\infty)$ must be at least $2$ .

If the conditions of a case are not satisﬁed then the case will be attempted as the algorithm guarantees that there will be no Liouvillian solution. However if the conditions are satisﬁed, this does not necessarily mean a solution exists. As an example $y^{″} = 1 / x^{6} y$ satisﬁes only case one, but running the algorithm on case one shows that there is no Liouvillian solution.

The following table summarizes the above conditions for each case.


Case	Allowed pole order for $r$	Allowed value for $O (\infty)$

1	${0, 1, 2, 4, 6, 8, \dots}$	${\dots, - 6, - 4, - 2, 0, 2, 3, 4, 5, 6, \dots}$

2	Need to have at least one pole of order $2$ or pole of odd order greater than $2$ . Any other pole order is allowed as long as the above condition is satisﬁed. The following are examples of pole orders which are allowed. ${1, 2}$ , ${1, 3}$ , ${2}$ , ${3}$ , ${3, 4}$ , ${1, 2, 5}$ .	no conditions

3	${1, 2}$	${2, 3, 4, 5, 6, 7, \dots}$

Table 1: Necessary conditions for each Kovacic case

Some observations: In case one, no odd order pole is allowed except for order 1. Case one is the only case that could have no pole in $r$ , which is the same as a pole of order zero. Case two and three require at least one pole. For case three, only poles of order $1$ or $2$ are allowed. If $O (\infty)$ is zero, then only possibility is either case one or two. For case one, if $O (\infty)$ is negative, then it must be even.

The above table also shows that when $r$ has only one pole of order $2$ and $O (\infty)$ equals $2$ or higher then all three cases are possible. Also, if $r$ has two poles one of order $1$ and the other of order $2$ and $O (\infty)$ equals $2$ or higher then all three cases are possible.

These are the only two possibilities where all three cases have the same necessary conditions.

2 Description of algorithm for each case

2.1 Case one

2.1.1 step 1

Assuming that the necessary conditions for case one are satisﬁed and $z^{″} = r z, r = \frac{s}{t}$ . Let $Γ$ be the set of all poles of $r$ . For each pole $c$ in this set, three quantities are calculated: Rational function ${[\sqrt{r}]}_{c}$ and two complex numbers $α_{c}^{+}, α_{c}^{-}$ .

How this is done depends on the order of the pole as described below. If the set $Γ$ is empty (when there are no poles), then this part is skipped.

If the pole $c$ has order $1$ then $\begin{aligned} {[\sqrt{r}]}_{c} & = 0 \\ α_{c}^{+} & = 1 \\ α_{c}^{-} & = 1 \end{aligned}$
If the pole $c$ is of order $2$ then $\begin{aligned} {[\sqrt{r}]}_{c} & = 0 \\ α_{c}^{+} & = \frac{1}{2} + \frac{1}{2} \sqrt{1 + 4 b} \\ α_{c}^{-} & = \frac{1}{2} - \frac{1}{2} \sqrt{1 + 4 b} \end{aligned}$

Where $b$ is the coeﬃcient of $\frac{1}{(x - c)^{2}}$ in the partial fraction decomposition of $r$ .
If the pole is of order ${4, 6, 8, \dots}$ (poles must be all even from the conditions of case one), then the computation is more involved. Let $2 v$ be the order of the pole. Hence if the pole was order 4, then $v = 2$ . Let ${[\sqrt{r}]}_{c}$ be the sum of terms involving $\frac{1}{(x - c)^{i}}$ for $2 \leq i \leq v$ in the Laurent series expansion of $\sqrt{r}$ (not $r$ ) at $c$ . Therefore $\begin{aligned} {[\sqrt{r}]}_{c} & = \sum_{i = 2}^{v} \frac{a_{i}}{{(x - c)}^{i}} \\ (1) & = \frac{a_{2}}{{(x - c)}^{2}} + \frac{a_{3}}{{(x - c)}^{3}} + \dots + \frac{a_{v}}{{(x - c)}^{v}} \end{aligned}$

$α_{c}^{+}, α_{c}^{-}$ are found using $\begin{aligned} α_{c}^{+} & = \frac{1}{2} (\frac{b}{a_{v}} + v) \\ α_{c}^{-} & = \frac{1}{2} (- \frac{b}{a_{v}} + v) \end{aligned}$

Where in the above $a_{v}$ is the coeﬃcient of the term $\frac{a_{v}}{(x - c)^{v}}$ in (1) and $b$ is the coeﬃcient of the term $\frac{1}{(x - c)^{v + 1}}$ in $r$ itself (found from the partial fraction decomposition), minus the coeﬃcient of same term in the Laurent series expansion of $\sqrt{r}$ at $c$ .

The coeﬃcients in the Laurent series can be obtained as follows. Given $r (x)$ with a pole of ﬁnite order $N$ at $x = c$ , then its Laurent series expansion at $c$ is given by the sum of the analytic part and the principal part of the of the Laurent series. The coeﬃcients $b_{n}$ are contained in the principal part of the series. $\begin{aligned} (2) & r (x) & = \sum_{n = 0}^{\infty} a_{n} {(x - c)}^{n} + \sum_{n = 1}^{N} \frac{b_{n}}{{(x - c)}^{n}} \\ = \sum_{n = 0}^{\infty} a_{n} {(x - c)}^{n} + \frac{b_{1}}{(x - c)} + \frac{b_{2}}{{(x - c)}^{2}} + \frac{b_{3}}{{(x - c)}^{3}} + \dots + \frac{b_{N}}{{(x - c)}^{N}} \end{aligned}$

To obtain $b_{1}$ (which is the residue of $r (x)$ at $c$ ), both sides of the above are multiplied by ${(x - c)}^{N}$ which gives $\begin{matrix} (3) & {(x - c)}^{N} r (x) = \sum_{n = 0}^{\infty} a_{n} {(x - c)}^{n + N} + b_{1} {(x - c)}^{N - 1} + b_{2} {(x - c)}^{N - 2} + \dots + b_{N} \end{matrix}$ Diﬀerentiating both sides of (3) $(N - 1)$ times w.r.t. $x$ gives $\frac{d^{N - 1}}{d x^{(N - 1)}} ({(x - c)}^{N} f (x)) = \sum_{n = 0}^{\infty} \frac{d^{N - 1}}{d x^{(N - 1)}} (a_{n} {(x - c)}^{n + N}) + b_{1} (N - 1)!$ Evaluating the above at $x = c$ gives $b_{1} = \frac{lim_{x \to c} \frac{d^{N - 1}}{d x^{(N - 1)}} ({(x - c)}^{N} r (x))}{(N - 1)!}$ To ﬁnd the next coeﬃcient $b_{2}$ , both sides of (3) are diﬀerentiated $(N - 2)$ times $\frac{d^{N - 2}}{d x^{(N - 2)}} ({(x - c)}^{N} r (x)) = \sum_{n = 0}^{\infty} \frac{d^{N - 2}}{d x^{(N - 2)}} (a_{n} {(x - c)}^{n + N}) + b_{1} (N - 1)! (x - c) + b_{2} (N - 2)!$ Evaluating the above at $x = c$ gives $b_{2} = \frac{lim_{x \to c} \frac{d^{N - 2}}{d x^{(N - 2)}} ({(x - c)}^{N} r (x))}{(N - 2)!}$ The above is repeated to ﬁnd $b_{3}, b_{4}, \dots, b_{N}$ . The general formula for ﬁnd coeﬃcient $b_{n}$ is therefore $\begin{array}{r} (4) & b_{n} = \frac{lim_{x \to c} \frac{d^{N - n}}{d x^{(N - n)}} ({(x - c)}^{N} r (x))}{(N - n)!} \end{array}$

For the special case of the last term $b_{N}$ the above simpliﬁes to $\begin{array}{r} (5) & b_{N} = lim_{x \to c} (x - c)^{N} r (x) \end{array}$

The above is implemented in the function laurent_coeff() in the Kovacic class.

This completes ﬁnding all the quantities ${{[\sqrt{r}]}_{c}, α_{c}^{+}, α_{c}^{+}}$ for each pole in the set $Γ$ for case one.

The next step calculates the following three quantities for $O (\infty)$ .

If $O (\infty) \leq 0$ , which must be even, then let $- 2 v = O (\infty)$ and ${[\sqrt{r}]}_{\infty}$ is then the sum of all terms $x^{i}$ for for $0 \leq i \leq v$ in the Laurent series expansion of $\sqrt{r}$ at $\infty$ . $\begin{aligned} (6) & {[\sqrt{r}]}_{\infty} & = \sum_{i = 0}^{v} a_{i} x^{v} = a_{0} + a_{1} x + a_{2} x^{2} \dots + a_{v} x^{v} \end{aligned}$

The coeﬃcients $a_{i}$ are found by setting $x = \frac{1}{y}$ in $r$ and then ﬁnding the Laurent series of $[\sqrt{r (y)}]$ expanded around $y = z e r o$ . The process for ﬁnding the coeﬃcient is the same one used as described earlier where now the limit is taken as $y$ approaches zero from the right. This gives all the terms of (6). This is implemented in the function laurent_coeff() in the Kovacic class.

The corresponding ${α_{\infty}^{+}, α_{\infty}^{-}}$ are given by $\begin{aligned} α_{\infty}^{+} & = \frac{1}{2} (\frac{b}{a_{v}} - v) \\ α_{\infty}^{-} & = \frac{1}{2} (- \frac{b}{a_{v}} - v) \end{aligned}$

Where $a_{v}$ is coeﬃcient of $x^{v}$ in (6) and $b$ is the coeﬃcient of $x^{v - 1}$ in $r$ itself (found using long division) minus the coeﬃcient of $x^{v - 1}$ in ${({[\sqrt{r}]}_{\infty})}^{2}$ .
If $O (\infty) = 2$ then ${[\sqrt{r}]}_{\infty} = 0$ . The corresponding ${α_{\infty}^{+}, α_{\infty}^{-}}$ are given by $\begin{aligned} α_{\infty}^{+} & = \frac{1}{2} + \frac{1}{2} \sqrt{1 + 4 b} \\ α_{\infty}^{-} & = \frac{1}{2} - \frac{1}{2} \sqrt{1 + 4 b} \end{aligned}$

Here $b = \frac{lcoef (s)}{lcoeff (t)}$ where $r = \frac{s}{t}$ . $lcoef (s)$ is the leading coeﬃcient of $s$ and similarly, $lcoef (t)$ is the leading coeﬃcient of $t$ .
If $O (\infty) > 2$ then $\begin{aligned} {[\sqrt{r}]}_{\infty} & = 0 \\ α_{\infty}^{+} & = 0 \\ α_{\infty}^{-} & = 1 \end{aligned}$

2.1.2 step 2

Using quantities calculated in step $1$ , the algorithm now searches for a non-negative integer $d$ using $\begin{aligned} d & = α_{\infty}^{\pm} - \sum_{c \in Γ} α_{c}^{\pm} \end{aligned}$

If non-negative $d$ is found, a candidate $ω_{d}$ is calculated using $\begin{aligned} ω_{d} & = \sum_{c \in Γ} ((\pm) {[\sqrt{r}]}_{c} + \frac{α_{c}^{\pm}}{x - c}) + (\pm) {[\sqrt{r}]}_{\infty} \end{aligned}$

If no non-negative integer $d$ could be found, then no Liouvillian solution exists using this case. Case two or three are tried next if these are available.

2.1.3 step 3

In this step the algorithm ﬁnds polynomial $p (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + x^{d}$ of degree $d$ . This is done by solving for the coeﬃcients $a_{i}$ from $\begin{aligned} (7) & p^{''} + 2 ω p^{'} + (ω^{'} + ω^{2} - r) p & = 0 \end{aligned}$

Where $ω$ is from the second step above and $r$ is from $z^{″} = r z$ .

For an example, if $d = 2$ , then $p (x) = x^{2} + a_{1} x + a_{0}$ is substituted in (3) and $a_{0}, a_{1}$ are solved for. If solution exists, then the solution to $z^{″} = r z$ will be $\begin{aligned} z & = p (x) e^{\int ω d x} \end{aligned}$

If the degree $d = 1$ then $p (x) = x + a_{0}$ and the same process is applied. If the degree $d = 0$ , then $p (x) = 1$ .

The ﬁrst basis solution to the original ode is now be found from $\begin{aligned} y_{1} & = z e^{- \frac{1}{2} \int a d x} \end{aligned}$

And the second basis solution using reduction of order formula is $\begin{aligned} y_{2} & = y_{1} \int \frac{e^{- \int a d x}}{y_{1}^{2}} d x \end{aligned}$

Hence the general solution to the original ode is $\begin{array}{r} y (x) = c_{1} y_{1} + c_{2} y_{2} \end{array}$

This completes the full algorithm for case 1. The part that needs most care is in ﬁnding ${{[\sqrt{r}]}_{c}, α_{c}^{\pm}, {[\sqrt{r}]}_{\infty}, α_{\infty}^{\pm}}$ . Once these are calculated, the rest of the algorithm is much more direct.

2.1.4 Algorithm ﬂow chart for case one

pict — Figure 1: Case 1 Kovacic algorithm

2.2 Case two

2.2.1 step 1

Assuming that the necessary conditions for case two are satisﬁed and $z^{″} = r z, r = \frac{s}{t}$ . Let $Γ$ be the set of all poles of $r$ . For each pole $c$ in this set, $E_{c}$ is found as follows

If the pole $c$ has order $1$ then $E_{c} = {4}$ .
If the pole $c$ is of order $2$ then $E_{c} = {2, 2 + 2 \sqrt{1 + 4 b}, 2 - 2 \sqrt{1 + 4 b}}$ where $b$ is the coeﬃcient of $\frac{1}{(x - c)^{2}}$ in the partial fraction decomposition of $r$ . In the above set $E_{c}$ , only integer values are kept.
If the pole $c$ is of order $v > 2$ then $E_{c} = {v}$

The next step is to determine $E_{\infty}$ .

If $O (\infty) > 2$ then $E_{\infty} = {0, 2, 4}$
If $O (\infty) = 2$ then $E_{\infty} = {2, 2 + 2 \sqrt{1 + 4 b}, 2 - 2 \sqrt{1 + 4 b}}$ where $b = \frac{lcoef (s)}{lcoeff (t)}$ where $r = \frac{s}{t}$ . $lcoef (s)$ is the leading coeﬃcient of $s$ and similarly $lcoef (t)$ is the leading coeﬃcient of $t$ . In the above set $E_{\infty}$ only integer values are kept.
If $O (\infty) < 2$ then $E_{\infty} = O (\infty)$ .

2.2.2 step 2

Using quantities calculated in step $1$ , the algorithm now searches for a non-negative integer $d$ using $\begin{aligned} d & = \frac{1}{2} (e_{\infty} - \sum_{c \in Γ} e_{c}) \end{aligned}$

Where in the above $e_{c} \in E_{c}$ , $e_{\infty} \in E_{\infty}$ found in step $1$ . If non-negative $d$ is found, then $\begin{aligned} θ & = \frac{1}{2} \sum_{c \in Γ} \frac{e_{c}}{x - c} \end{aligned}$

If no non-negative integer $d$ could be found, then no Liouvillian solution exists using this case. Case three is tried next if it is available.

2.2.3 step 3

In this step the algorithm determines a polynomial $p (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + x^{d}$ of degree $d$ . This is done by solving for the coeﬃcients $a_{i}$ from $\begin{array}{r} (1) & p^{‴} + 3 θ p^{″} + (3 θ^{2} + 3 θ^{'} - 4 r) p^{'} + (θ^{″} + 3 θ θ^{'} + θ^{3} - 4 r θ - 2 r^{'}) p = 0 \end{array}$

Where $θ$ was found in step $2$ and $r$ is from $z^{″} = r z$ . If $p (x)$ can be found that satisﬁes (1) then $\begin{array}{r} (2) & ϕ = θ + \frac{p^{'}}{p} \end{array}$

$ω$ is then solved for from $\begin{aligned} (3) & ω^{2} - ϕ ω + (\frac{1}{2} ϕ^{'} + \frac{1}{2} ϕ^{2} - r) & = 0 \end{aligned}$

If solution $ω$ to (3) can be found, then the solution to $z^{″} = r z$ is given by $\begin{array}{r} z = e^{\int ω d x} \end{array}$

This completes the full algorithm for case two. The general solution to the original ode is now determined as outlined at the end of case one above.

2.2.4 Algorithm ﬂow chart for case two

2.3 Case three

2.3.1 step 1

Assuming the necessary conditions for case three are satisﬁed and $z^{″} = r z, r = \frac{s}{t}$ . Let $Γ$ be the set of all poles of $r$ . Recall that case three can have either a pole of order 1 or order 2 only. For each pole $c$ in this set, $E_{c}$ is found as follows

If the pole $c$ has order $1$ then $E_{c} = {12}$ .
If the pole $c$ is of order $2$ then $\begin{aligned} (1) & E_{c} & = {6 + \frac{12 k}{n} \sqrt{1 + 4 b}} for k & = - \frac{n}{2} \dots \frac{n}{2} \end{aligned}$

Where $k$ is incremented by $1$ each time, and $n$ is any of ${4, 6, 12}$ and $b$ is the coeﬃcient of $\frac{1}{(x - c)^{2}}$ in the partial fraction decomposition of $r$ . In the above set $E_{c}$ , only integer values are kept. For an example, when $n = 4$ then $k = {- 2, - 1, 0, 1, 2}$ and $E_{c} = {6 - 6 \sqrt{1 + 4 b}, 6 - 3 \sqrt{1 + 4 b}, 6, 6 + 3 \sqrt{1 + 4 b}, 6 + 6 \sqrt{1 + 4 b}}$ and similarly for $n = 6$ and $n = 12$ .

The next step determines $E_{\infty}$ . This is found using same formula as (1) but $b$ is calculated diﬀerently using $b = \frac{lcoef (s)}{lcoeff (t)}$ where $r = \frac{s}{t}$ . $lcoef (s)$ is the leading coeﬃcient of $s$ and $lcoef (t)$ is the leading coeﬃcient of $t$ .

2.3.2 step 2

Using quantities calculated in step $1$ , the algorithm now searches for a non-negative integer $d$ using $\begin{aligned} d & = \frac{n}{12} (e_{\infty} - \sum_{c \in Γ} e_{c}) \end{aligned}$

Where in the above $e_{c} \in E_{c}$ , $e_{\infty} \in E_{\infty}$ $n$ is any of ${4, 6, 12}$ values. If non-negative $d$ is found, then $\begin{aligned} θ & = \frac{n}{12} \sum_{c \in Γ} \frac{e_{c}}{x - c} \end{aligned}$

The sum above is over all families of ${e_{\infty}, e_{c}}$ which generated the non-negative integer $d$ . Next deﬁne $\begin{aligned} S & = \prod_{c \in Γ} (x - c) \end{aligned}$

The product above is over families of ${e_{\infty}, e_{c}}$ which generated the non-negative integer $d$ . If no non-negative integer $d$ is found, then no Liouvillian solution exists.

2.3.3 step 3

In this step the algorithm determines a polynomial $p (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + x^{d}$ of degree $d$ . Deﬁne set of polynomials ${P_{n}, P_{n - 1}, \dots, P_{- 1}$ where $\begin{aligned} P_{n} & = - p (x) \\ P_{i - 1} & = - S P_{i}^{'} + ((n - i) S^{'} - S θ) P_{i} - (n - i) (i + 1) S^{2} r P_{i + 1} i = n \dots 0 \end{aligned}$

The last polynomial $P_{- 1} (x)$ is used to solve for the coeﬃcients $a_{i}$ using $\begin{aligned} (2) & P_{- 1} (x) & = 0 \end{aligned}$

In Maple this is done using the solve command with the identity option. If it is possible to ﬁnd coeﬃcients $a_{i}$ such that (2) is satisﬁed, then deﬁne the equation $\begin{array}{r} \sum_{i = 0}^{n} \frac{S^{i} P_{i} (x)}{(n - i)!} ω^{i} = 0 \end{array}$

$ω$ is solved for from the above equation. If solution $ω$ is found then the solution to $z^{″} = r z$ will be $\begin{aligned} z & = e^{\int ω d x} \end{aligned}$

This completes the full algorithm for case three. The general solution to the original ode can now be determined as outlined at the end of case one above.

2.3.4 Algorithm ﬂow chart for case three

2.4 Statistics and discussion of results obtained using Kovacic algorithm

This gives summary of results obtained using testsuite of $3000$ diﬀerential equations, all of which were selected as linear with rational coeﬃcients as functions of $x$ that can be solved using this algorithm.

The ode’s used in the testsuite were collected by the author and stored in sql database. These were collected from a number of standard textbooks and other references such as “Diﬀerential Equations. E. Kamke. 3th edition. Chelsea.” and “Ordinary Diﬀerential Equations And Their Solutions. Murphy, George Moseley. Dover. 2011”.

All the ode’s were successfully solved using the Kovacic algorithm as implemented here and each solution was veriﬁed using Maple odetest.

The following diagram shows the percentage of ode’s solved using each case.

Case $3$ was required for solving only $3$ odes. It used $n = 4$ for all $3$ ode’s. $n = 6$ and $n = 12$ were not reached or required to try. Recall that $n$ for case $3$ is the degree of the polynomial in $ω$ used to solve for in order to ﬁnd the $z$ solution from $z = e^{\int ω d x}$ .

This result shows that case $1$ and $2$ combined is all what is needed to solve $99.9$ % of ode’s used in practice. Larger collection of ode’s than the $3000$ used could produce diﬀerent results, but the overall trend is that case $3$ is rarely needed in practice and within case $3$ , $n = 6$ and $n = 12$ are even less likely to be required.

When forcing the algorithm to use case $3$ and only use $n = 12$ , this resulted in a very long computation time on some ode’s. For an example, using ode $y^{″} + x y^{'} + y = 0$ which satisﬁes all three cases, and asking the solver to use case $3$ and $n = 12$ , it was found that it required $p (x)$ of degree $d = 24$ in order to ﬁnd $ω$ of degree $12$ that can be solved. The total number of trials in step 3 of case three to ﬁnd such solution was found to be $2367$ . This took over 30 minutes to complete.

In comparison, the same ode was solved using case one in less than one second giving the same solution on the same computer.

The testsuite also calculates the distribution of cases which has its necessary conditions satisﬁed for each ode. Recall that having the necessary conditions for a case satisﬁed does not mean a solution would be found using that case. The following bar chart shows the percentages of the $3000$ ode’s that satisﬁed the necessary conditions each case. This chart shows that many ode’s satisfy the conditions for more than one case at the same time.

3 Worked example for each case

3.1 case one

3.1.1 Example 1

Given the ode $\begin{array}{r} (2 x + 1) y^{″} - 2 y^{'} - (2 x + 3) y = 0 \end{array}$

Converting it $y^{″} + a y^{'} + b y = 0$ gives $\begin{array}{r} y^{″} - \frac{2}{2 x + 1} y^{'} - \frac{2 x + 3}{2 x + 1} y = 0 \end{array}$

Where $a = - \frac{2}{2 x + 1}, b = - \frac{2 x + 3}{2 x + 1}$ . Applying the transformation $z = y e^{\frac{1}{2} \int a d x}$ gives $z^{″} = r z$ where $r = \frac{1}{4} a^{2} + \frac{1}{2} a^{'} - b$ . This results in $\begin{aligned} r & = \frac{s}{t} \\ = \frac{4 x^{2} + 8 x + 6}{(2 x + 1)^{2}} = \frac{s}{t} \end{aligned}$

There is one pole at $x = - \frac{1}{2}$ , hence $Γ = {- \frac{1}{2}}$ . The order is $2$ and $O (\infty) = \deg (t) - \deg (s) = 0$ . Table 1 shows that the necessary conditions for case one and two are both satisﬁed. This is solved ﬁrst using case one. Since the order of the pole is $2$ , then $\begin{aligned} (1) & {[\sqrt{r}]}_{c} & = 0 \\ α_{c}^{+} & = \frac{1}{2} + \frac{1}{2} \sqrt{1 + 4 b} \\ α_{c}^{-} & = \frac{1}{2} - \frac{1}{2} \sqrt{1 + 4 b} \end{aligned}$

Where $b$ is the coeﬃcient of $\frac{1}{(x - c)^{2}} = \frac{1}{(x + \frac{1}{2})^{2}}$ in the partial fraction decomposition of $r$ which is $r = 1 + \frac{3}{4} \frac{1}{(x + \frac{1}{2})^{2}} + \frac{1}{x + \frac{1}{2}}$ . Therefore $b = \frac{3}{4}$ and the above becomes $\begin{aligned} {[\sqrt{r}]}_{c} & = 0 \\ α_{c}^{+} & = \frac{1}{2} + \frac{1}{2} \sqrt{1 + 4 (\frac{3}{4})} = \frac{3}{2} \\ α_{c}^{-} & = \frac{1}{2} - \frac{1}{2} \sqrt{1 + 4 (\frac{3}{4})} = - \frac{1}{2} \end{aligned}$

Since $O (\infty) = 0$ then $v = 0$ and ${[\sqrt{r}]}_{\infty}$ is the sum of all terms $x^{i}$ for $0 \leq i \leq v$ in the Laurent series expansion of $\sqrt{r}$ at $\infty$ which is found as follows. Since $\sqrt{r} = \sqrt{\frac{4 x^{2} + 8 x + 6}{(2 x + 1)^{2}}}$ then setting $x = \frac{1}{y}$ gives $\sqrt{r (y)} = \sqrt{\frac{4 {(\frac{1}{y})}^{2} + 8 \frac{1}{y} + 6}{{(2 \frac{1}{y} + 1)}^{2}}}$ and since $v = 0$ then the constant term is $lim_{y \to 0} \sqrt{r (y)} = 1$ . Therefore ${[\sqrt{r (x)}]}_{\infty} = 1$ . Hence $a = 1$ .

$b$ is the coeﬃcient of $x^{v - 1} = \frac{1}{x}$ in $r$ minus the coeﬃcient of $\frac{1}{x}$ in ${({[\sqrt{r}]}_{\infty})}^{2} = 1$ which is zero since there is no term $\frac{1}{x}$ . Because $v = 0$ , long division is used is used ﬁnd the coeﬃcient of $\frac{1}{x}$ in $r$ . $\begin{aligned} r & = \frac{s}{t} \\ = \frac{4 x^{2} + 8 x + 6}{4 x^{2} + 4 x + 1} \\ = Q + \frac{R}{t} \\ = 1 + \frac{4 x + 5}{4 x^{2} + 4 x + 1} \end{aligned}$

The coeﬃcient of $\frac{1}{x}$ in $r$ is the leading coeﬃcient in $R$ minus the leading coeﬃcient in $t$ which gives $\frac{4}{4} = 1$ . Therefore $b = 1 - 0 = 1$ . This results in $\begin{aligned} {[\sqrt{r}]}_{\infty} & = 1 \\ α_{\infty}^{+} & = \frac{1}{2} (\frac{b}{a_{v}} - v) & = \frac{1}{2} (\frac{1}{1} - 0) & = \frac{1}{2} \\ α_{\infty}^{-} & = \frac{1}{2} (- \frac{b}{a_{v}} - v) & = \frac{1}{2} (- \frac{1}{1} - 0) & = - \frac{1}{2} \end{aligned}$

The above completes step 1 for case one. Step $2$ searches for a non-negative integer $d$ using $\begin{aligned} (2) & d & = α_{\infty}^{\pm} - \sum_{c \in Γ} α_{c}^{\pm} \end{aligned}$

Where the above is carried over all possible combinations resulting in the following 4 possibilities (in this example, there is only one pole, hence the sum contains only one term) and the number of possible combinations is therefore $2^{2} = 4$ $\begin{aligned} d & = α_{\infty}^{+} - α_{c}^{+} & = \frac{1}{2} - (\frac{3}{2}) & = - 1 \\ d & = α_{\infty}^{+} - α_{c}^{-} & = \frac{1}{2} - (- \frac{1}{2}) & = 1 \\ d & = α_{\infty}^{-} - α_{c}^{+} & = - \frac{1}{2} - (\frac{3}{2}) & = - 2 \\ d & = α_{\infty}^{-} - α_{c}^{-} & = - \frac{1}{2} - (- \frac{1}{2}) & = 0 \end{aligned}$

The above shows there are two possible $d$ values to use. $d = 1$ or $d = 0$ . Each is tried until one produces a solution or all fail to do so. For each valid $d$ found an $ω$ is found using $\begin{aligned} ω & = \sum_{c} ((\pm) {[\sqrt{r}]}_{c} + \frac{α_{c}^{\pm}}{x - c}) + (\pm) {[\sqrt{r}]}_{\infty} \end{aligned}$

But ${[\sqrt{r}]}_{c} = 0$ in this example, hence the above simpliﬁes to $\begin{aligned} ω & = \sum_{c} (\frac{α_{c}^{\pm}}{x - c}) + (\pm) {[\sqrt{r}]}_{\infty} \end{aligned}$

Since there is one pole, then the candidate $ω$ to try are the following $\begin{aligned} ω & = \frac{α_{c}^{+}}{x - c} + (+ 1) {[\sqrt{r}]}_{\infty} \\ ω & = \frac{α_{c}^{+}}{x - c} + (- 1) {[\sqrt{r}]}_{\infty} \\ ω & = \frac{α_{c}^{-}}{x - c} + (+ 1) {[\sqrt{r}]}_{\infty} \\ ω & = \frac{α_{c}^{-}}{x - c} + (- 1) {[\sqrt{r}]}_{\infty} \end{aligned}$

Substituting the known values found in step 1 into the above gives $\begin{aligned} ω & = \frac{\frac{3}{2}}{x + \frac{1}{2}} + (+ 1) (1) = \frac{6 + 2 x}{2 x + 3} \\ ω & = \frac{\frac{3}{2}}{x + \frac{1}{2}} + (- 1) (1) = - \frac{2 x}{2 x + 3} \\ ω & = \frac{- \frac{1}{2}}{x + \frac{1}{2}} + (+ 1) (1) = \frac{2 x}{2 x + 1} \\ ω & = \frac{- \frac{1}{2}}{x + \frac{1}{2}} + (- 1) (1) = - \frac{2 (1 + x)}{2 x + 1} \end{aligned}$

So there are two possible $d$ values to try, and for each, there are $4$ possible $w (x)$ , which gives $8$ possible tries. This completes step 2. For each trial, step 3 is now invoked.

Starting with $d = 0$ and using $ω = \frac{2 x}{2 x + 1}$ , and since the degree is $d = 0$ then $p (x) = 1$ . This polynomial is now checked to see if it satisﬁes $\begin{aligned} (3) & p^{″} + 2 ω p^{'} + (ω^{'} + ω^{2} - r) p & = 0 \\ (ω^{'} + ω^{2} - r) p & = 0 \\ \frac{d}{d x} (\frac{2 x}{2 x + 1}) + {(\frac{2 x}{2 x + 1})}^{2} - \frac{4 x^{2} + 8 x + 6}{(2 x + 1)^{2}} & = 0 \\ - \frac{4}{2 x + 1} & = 0 \end{aligned}$

Since the left side is not identically zero, then this candidate $ω$ has failed. Carrying out this process for the other 3 possible $ω$ values shows that non are satisﬁed as well. $d = 1$ is now tried. This implies the polynomial is $p (x) = a_{0} + x$ . The coeﬃcient $a_{0}$ needs to be determined such that $p^{″} + 2 ω p^{'} + (ω^{'} + ω^{2} - r) p = 0$ is satisﬁed. Starting with $ω = \frac{2 x}{2 x + 1}$ gives $\begin{aligned} p^{″} + 2 ω p^{'} + (ω^{'} + ω^{2} - r) p & = 0 \\ 2 ω p^{'} + (ω^{'} + ω^{2} - r) p & = 0 \end{aligned}$

Substituting $p = a_{0} + x$ and $ω = \frac{2 x}{2 x + 1}$ and $r = \frac{4 x^{2} + 8 x + 6}{(2 x + 1)^{2}}$ into the above and simplifying gives $- \frac{4 a_{0}}{2 x + 1} = 0$ . This implies that (3) can be satisﬁed for $a_{0} = 0$ . Therefore the polynomial of degree one is found and given by $\begin{aligned} p (x) & = x \end{aligned}$

Therefore the solution to $z^{″} = r z$ is $\begin{aligned} z & = p e^{\int ω d x} \\ = x e^{\int \frac{2 x}{2 x + 1} d x} \\ = x e^{x - \frac{\ln (2 x + 1)}{2}} \\ = \frac{x e^{x}}{\sqrt{2 x + 1}} \end{aligned}$

Given this solution for $z (x)$ , the ﬁrst basis solution of the original ode in $y$ is found using the inverse of the original transformation used to generate the $z$ ode which is $z = y e^{\frac{1}{2} \int a d x}$ , therefore $\begin{aligned} y_{1} & = z e^{- \frac{1}{2} \int a d x} \\ = \frac{x e^{x}}{\sqrt{2 x + 1}} e^{- \frac{1}{2} \int - \frac{2}{2 x + 1} d x} \\ = \frac{x e^{x}}{\sqrt{2 x + 1}} e^{\frac{\ln (2 x + 1}{2}} \\ = \frac{x e^{x}}{\sqrt{2 x + 1}} \sqrt{2 x + 1} \\ = x e^{x} \end{aligned}$

The second basis solution is found using reduction of order $\begin{aligned} y_{2} & = y_{1} \int \frac{e^{\int - a d x} d x}{y_{1}^{2}} d x \\ = x e^{x} \int \frac{e^{\int - \frac{- 2}{2 x + 1} d x} d x}{(x e^{x})^{2}} d x \\ = x e^{x} \int \frac{e^{\ln (2 x + 1)}}{(x e^{x})^{2}} d x \\ = - e^{- x} \end{aligned}$

Therefore the general solution to the original ode $(2 x + 1) y^{″} - 2 y^{'} - (2 x + 3) y = 0$ is $\begin{aligned} y (x) & = c_{1} y_{1} + x_{2} y_{2} \\ y (x) & = c_{1} x e^{x} - c_{2} e^{- x} \end{aligned}$

This completes the solution.

3.1.2 Example 2

Given the ode $\begin{array}{r} x^{2} (x^{2} - 2 x + 1) y^{″} - x (3 + x) y^{'} + (4 + x) y = 0 \end{array}$

Converting it $y^{″} + a y^{'} + b y = 0$ gives $\begin{array}{r} y^{″} - \frac{x (3 + x)}{x^{2} (x^{2} - 2 x + 1)} y^{'} + \frac{4 + x}{x^{2} (x^{2} - 2 x + 1)} y = 0 \end{array}$

Where $a = - \frac{x (3 + x)}{x^{2} (x^{2} - 2 x + 1)}, b = \frac{4 + x}{x^{2} (x^{2} - 2 x + 1)}$ . Applying the transformation $z = y e^{\frac{1}{2} \int a d x}$ gives $z^{″} = r z$ where $r = \frac{1}{4} a^{2} + \frac{1}{2} a^{'} - b$ . This results in $\begin{aligned} r & = \frac{s}{t} \\ = \frac{7 x^{2} + 10 x - 1}{4 x^{2} {(x - 1)}^{4}} \end{aligned}$

There is one pole at $x = 0$ of order $2$ and pole at $x = 1$ of order $4$ , hence $Γ = {0, 1}$ , and $O (\infty) = \deg (t) - \deg (s) = 6 - 2 = 4$ . Table 1 shows that the necessary conditions for case one and two are both satisﬁed. For the pole at $0$ and since its order is $2$ then $\begin{aligned} (1) & {[\sqrt{r}]}_{0} & = 0 \\ α_{0}^{+} & = \frac{1}{2} + \frac{1}{2} \sqrt{1 + 4 b} \\ α_{0}^{-} & = \frac{1}{2} - \frac{1}{2} \sqrt{1 + 4 b} \end{aligned}$

Where $b$ is the coeﬃcient of $\frac{1}{(x - 0)^{2}} = \frac{1}{x^{2}}$ in the partial fraction decomposition of $r$ which is $\begin{aligned} (2) & r & = \frac{4}{{(x - 1)}^{4}} - \frac{2}{{(x - 1)}^{3}} - \frac{1}{4 x^{2}} - \frac{3}{2 (x - 1)} + \frac{3}{2 x} + \frac{7}{4 {(x - 1)}^{2}} \end{aligned}$

The above shows that $b = - \frac{1}{4}$ . Equation (1) becomes becomes $\begin{aligned} {[\sqrt{r}]}_{0} & = 0 \\ α_{0}^{+} & = \frac{1}{2} + \frac{1}{2} \sqrt{1 - 4 \frac{1}{4}} = \frac{1}{2} \\ α_{0}^{-} & = \frac{1}{2} - \frac{1}{2} \sqrt{1 - 4 \frac{1}{4}} = \frac{1}{2} \end{aligned}$

For the second pole at $x = 1$ , since its order is $4$ , then $2 v = 4$ or $v = 2$ . Therefore the corresponding ${[\sqrt{r}]}_{1}$ is the sum of all terms involving $\frac{1}{(x - 1)^{i}}$ for $2 \leq i \leq v$ or $2 \leq i \leq 2$ in the Laurent series expansion of $\sqrt{r}$ (not $r$ ) around this pole. This results in $\begin{aligned} [\sqrt{r}]_{1} & = \sum_{i = 2}^{2} \frac{a_{i}}{(x - 1)^{i}} \\ (3) & = \frac{a_{2}}{(x - 1)^{2}} \end{aligned}$

$a_{2}$ is found using $\begin{aligned} a_{2} & = lim_{x \to 1} (x - 1)^{2} \sqrt{r} \\ = lim_{x \to 1} (x - 1)^{2} \sqrt{\frac{7 x^{2} + 10 x - 1}{4 x^{2} {(x - 1)}^{4}}} \\ = 2 \end{aligned}$

Therefore $\begin{aligned} [\sqrt{r}]_{1} & = \frac{2}{(x - 1)^{2}} \\ α_{1}^{+} & = \frac{1}{2} (\frac{b}{a} + v) & = \frac{1}{2} (\frac{b}{2} + 2) \\ (4) & α_{1}^{-} & = \frac{1}{2} (- \frac{b}{a} + v) & = \frac{1}{2} (- \frac{b}{2} + 2) \end{aligned}$

What remains is to determine $b$ . This is the coeﬃcient of $\frac{1}{(x - c)^{v + 1}} = \frac{1}{(x - 1)^{3}}$ in the partial fraction decomposition of $r$ which from (2) is $- 2$ minus the coeﬃcient of same term in $[\sqrt{r}]_{1}$ which from (3) is zero. Therefore $b = - 2 - 0 = - 2$ . (4) now becomes $\begin{aligned} [\sqrt{r}]_{1} & = \frac{2}{(x - 1)^{2}} \\ α_{1}^{+} & = \frac{1}{2} (\frac{b}{a} + v) & = \frac{1}{2} (\frac{- 2}{2} + 2) & = \frac{1}{2} \\ (5) & α_{1}^{-} & = \frac{1}{2} (- \frac{b}{a} + v) & = \frac{1}{2} (- \frac{- 2}{2} + 2) & = \frac{3}{2} \end{aligned}$

The above completes ﬁnding $[\sqrt{r}]_{c}, α_{c}^{+}, α_{c}^{+}$ for all poles in the set $Γ$ .

Since the order of $r$ at $\infty$ is $4 > 2$ then $\begin{aligned} [\sqrt{r}]_{\infty} & = 0 \\ α_{\infty}^{+} & = 0 \\ α_{\infty}^{-} & = 1 \end{aligned}$

This completes the ﬁrst step of the solution. The following tables summarizes the ﬁndings so far


pole $c$ location	pole order	$[\sqrt{r}]_{c}$	$α_{c}^{+}$	$α_{c}^{-}$

$0$	$2$	$0$	$\frac{1}{2}$	$\frac{1}{2}$

$1$	$4$	$\frac{2}{{(x - 1)}^{2}}$	$\frac{1}{2}$	$\frac{3}{2}$

Table 2: First step, case one.

Γ

set information


Order of $r$ at $\infty$	$[\sqrt{r}]_{\infty}$	$α_{\infty}^{+}$	$α_{\infty}^{-}$

$4$	$0$	$0$	$1$

Table 3: First step, case one.

O (\infty)

information

This completes step 1 for case one. Step 2 searches for non-negative integer $d$ using $\begin{aligned} d & = α_{\infty}^{\pm} - \sum_{c \in Γ} α_{c}^{\pm} \end{aligned}$

Where the above is carried over all possible combinations resulting in the following 8 possibilities (in this example, there are two poles, hence the sum contains two term) and the number of possible combinations is therefore $2^{3} = 8$ $\begin{aligned} d & = α_{\infty}^{+} - (α_{0}^{+} + α_{1}^{+}) & = 0 - (\frac{1}{2} + \frac{1}{2}) & = - 1 \\ d & = α_{\infty}^{+} - (α_{0}^{+} + α_{1}^{-}) & = 0 - (\frac{1}{2} + \frac{3}{2}) & = - 2 \\ d & = α_{\infty}^{+} - (α_{0}^{-} + α_{1}^{+}) & = 0 - (\frac{1}{2} + \frac{1}{2}) & = - 1 \\ d & = α_{\infty}^{+} - (α_{0}^{-} + α_{1}^{-}) & = 0 - (\frac{1}{2} + \frac{3}{2}) & = - 2 \\ d & = α_{\infty}^{-} - (α_{0}^{+} + α_{1}^{+}) & = 1 - (\frac{1}{2} + \frac{1}{2}) & = 0 \\ d & = α_{\infty}^{-} - (α_{0}^{+} + α_{1}^{-}) & = 1 - (\frac{1}{2} + \frac{3}{2}) & = - 1 \\ d & = α_{\infty}^{-} - (α_{0}^{-} + α_{1}^{+}) & = 1 - (\frac{1}{2} + \frac{1}{2}) & = 0 \\ d & = α_{\infty}^{-} - (α_{0}^{-} + α_{1}^{-}) & = 1 - (\frac{1}{2} + \frac{3}{2}) & = - 1 \end{aligned}$

There is only one possible $d = 0$ values to use. Candidate $ω$ are now found using $\begin{aligned} ω & = \sum_{c} ((\pm) {[\sqrt{r}]}_{c} + \frac{α_{c}^{\pm}}{x - c}) + (\pm) {[\sqrt{r}]}_{\infty} \end{aligned}$

Which gives $\begin{aligned} ω_{1} & = ((+) {[\sqrt{r}]}_{0} + \frac{α_{0}^{+}}{x}) + ((+) {[\sqrt{r}]}_{1} + \frac{α_{1}^{+}}{x - 1}) + (+) {[\sqrt{r}]}_{\infty} \\ ω_{2} & = ((+) {[\sqrt{r}]}_{0} + \frac{α_{0}^{+}}{x}) + ((-) {[\sqrt{r}]}_{1} + \frac{α_{1}^{-}}{x - 1}) + (+) {[\sqrt{r}]}_{\infty} \\ ω_{3} & = ((-) {[\sqrt{r}]}_{0} + \frac{α_{0}^{-}}{x}) + ((+) {[\sqrt{r}]}_{1} + \frac{α_{1}^{+}}{x - 1}) + (+) {[\sqrt{r}]}_{\infty} \\ ω_{4} & = ((-) {[\sqrt{r}]}_{0} + \frac{α_{0}^{-}}{x}) + ((-) {[\sqrt{r}]}_{1} + \frac{α_{1}^{-}}{x - 1}) + (+) {[\sqrt{r}]}_{\infty} \\ ω_{5} & = ((+) {[\sqrt{r}]}_{0} + \frac{α_{0}^{+}}{x}) + ((+) {[\sqrt{r}]}_{1} + \frac{α_{1}^{+}}{x - 1}) + (-) {[\sqrt{r}]}_{\infty} \\ ω_{6} & = ((+) {[\sqrt{r}]}_{0} + \frac{α_{0}^{+}}{x}) + ((-) {[\sqrt{r}]}_{1} + \frac{α_{1}^{-}}{x - 1}) + (-) {[\sqrt{r}]}_{\infty} \\ ω_{7} & = ((-) {[\sqrt{r}]}_{0} + \frac{α_{0}^{-}}{x}) + ((+) {[\sqrt{r}]}_{1} + \frac{α_{1}^{+}}{x - 1}) + (-) {[\sqrt{r}]}_{\infty} \\ ω_{8} & = ((-) {[\sqrt{r}]}_{0} + \frac{α_{0}^{-}}{x}) + ((-) {[\sqrt{r}]}_{1} + \frac{α_{1}^{-}}{x - 1}) + (-) {[\sqrt{r}]}_{\infty} \end{aligned}$

Substituting values found in step 1 into the above gives $\begin{aligned} ω_{1} & = (\frac{\frac{1}{2}}{x}) + ((+) \frac{2}{(x - 1)^{2}} + \frac{\frac{1}{2}}{x - 1}) & = \frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}} \\ ω_{2} & = (\frac{\frac{1}{2}}{x}) + ((-) \frac{2}{(x - 1)^{2}} + \frac{\frac{3}{2}}{x - 1}) & = \frac{4 x^{2} - 9 x + 1}{2 x {(x - 1)}^{2}} \\ ω_{3} & = (\frac{\frac{1}{2}}{x}) + ((+) \frac{2}{(x - 1)^{2}} + \frac{\frac{1}{2}}{x - 1}) & = \frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}} \\ ω_{4} & = (\frac{\frac{1}{2}}{x}) + ((-) \frac{2}{(x - 1)^{2}} + \frac{\frac{3}{2}}{x - 1}) & = \frac{4 x^{2} - 9 x + 1}{2 x {(x - 1)}^{2}} \\ ω_{5} & = (\frac{\frac{1}{2}}{x}) + ((+) \frac{2}{(x - 1)^{2}} + \frac{\frac{1}{2}}{x - 1}) & = \frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}} \\ ω_{6} & = (\frac{\frac{1}{2}}{x}) + ((-) \frac{2}{(x - 1)^{2}} + \frac{\frac{3}{2}}{x - 1}) & = \frac{4 x^{2} - 9 x + 1}{2 x {(x - 1)}^{2}} \\ ω_{7} & = (\frac{\frac{1}{2}}{x}) + ((+) \frac{2}{(x - 1)^{2}} + \frac{\frac{1}{2}}{x - 1}) & = \frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}} \\ ω_{8} & = (\frac{\frac{1}{2}}{x}) + ((-) \frac{2}{(x - 1)^{2}} + \frac{\frac{3}{2}}{x - 1}) & = \frac{4 x^{2} - 9 x + 1}{2 x {(x - 1)}^{2}} \end{aligned}$

Which shows there are only two diﬀerent $ω$ to try, these are $ω_{1}, ω_{2}$ . This complete step 2. For each trial, step 3 is now invoked. Starting with $d = 0$ and $ω = ω_{1} = \frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}}$ . Since the degree $d = 0$ then $p (x) = 1$ . This polynomial needs to satisfy $\begin{aligned} (6) & p^{″} + 2 ω p^{'} + (ω^{'} + ω^{2} - r) p & = 0 \\ (ω^{'} + ω^{2} - r) p & = 0 \\ \frac{d}{d x} (\frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}}) + {(\frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}})}^{2} - \frac{7 x^{2} + 10 x - 1}{4 x^{2} {(x - 1)}^{4}} & = 0 \\ 0 & = 0 \end{aligned}$

Because the equation is satisﬁed, the polynomial $p (x) = 1$ can be used. The solution to $z^{″} = r z$ is now found from $\begin{aligned} z & = p e^{\int ω d x} \\ = e^{\int \frac{2 x^{2} + x + 1}{2 x {(x - 1)}^{2}} d x} \\ = e^{\frac{\ln (x - 1)}{2} - \frac{2}{x - 1} + \frac{\ln (x)}{2}} \\ = \sqrt{x - 1} \sqrt{x} e^{- \frac{2}{x - 1}} \end{aligned}$

Given this solution for $z (x)$ , the ﬁrst basis solution of the original ode in $y$ is found using the inverse of the original transformation used to generate the $z$ ode which is $z = y e^{\frac{1}{2} \int a d x}$ , therefore $\begin{aligned} y_{1} & = z e^{- \frac{1}{2} \int a d x} \\ = \sqrt{x - 1} \sqrt{x} e^{- \frac{2}{x - 1}} e^{- \frac{1}{2} \int - \frac{x (3 + x)}{x^{2} (x^{2} - 2 x + 1)} d x} \end{aligned}$

Which simpliﬁes to $\begin{aligned} y_{1} & = \frac{x^{2}}{x - 1} e^{- \frac{4}{x - 1}} \end{aligned}$

The second solution $y_{2}$ to the original ode is found using reduction of order as was done in the ﬁrst example.

3.1.3 Example 3

This ode is a standard second order representing the oscillating harmonics ode with constant coeﬃcients and does not require Kovacic algorithm to solve it as it can be readily solved using standard method by ﬁnding the roots of the characteristic equation. It is included here in order to illustrate the Kovacic algorithm. $\begin{aligned} y^{″} + y^{'} + y & = 0 \\ A y^{″} + B y^{'} + C y & = 0 \end{aligned}$

Converting it to $z^{″} = r z$ as shown before gives $\begin{aligned} z^{″} & = \frac{s}{t} z \\ = \frac{- 3}{4} z \end{aligned}$

Hence $r = \frac{- 3}{4}$ . There are no poles therefore $Γ = {}$ , and $O (\infty) = \deg (t) - \deg (s) = 0$ . Table 1 shows that the necessary conditions for case one are only satisﬁed. Since the set $Γ$ is empty, then only the quantities related to $O (\infty)$ need to be calculated. The order of $r$ at $\infty$ is $O_{r} (\infty) = 0$ therefore $v = 0$ . $r$ has no $x$ in it, hence the Laurent series of $\sqrt{r}$ at $\infty$ is itself $\begin{array}{r} \sqrt{r} = \frac{i \sqrt{3}}{2} \end{array}$

Therefore $\begin{array}{r} a = \frac{i \sqrt{3}}{2} \end{array}$

And since $r$ is constant then $b = 0$ . Hence $\begin{aligned} [\sqrt{r}]_{\infty} & = \frac{i \sqrt{3}}{2} \\ α_{\infty}^{+} & = \frac{1}{2} (\frac{b}{a} - v) & = \frac{1}{2} (\frac{0}{\frac{i \sqrt{3}}{2}} - 0) & = 0 \\ α_{\infty}^{-} & = \frac{1}{2} (- \frac{b}{a} - v) & = \frac{1}{2} (- \frac{0}{\frac{i \sqrt{3}}{2}} - 0) & = 0 \end{aligned}$

This completes step 1 for case one. Step 2 searches for non-negative integer $d$ using $\begin{aligned} d & = α_{\infty}^{\pm} - \sum_{c \in Γ} α_{c}^{\pm} \end{aligned}$

Since there are no poles then $\begin{aligned} d & = α_{\infty}^{-} \\ = 0 \end{aligned}$

Since $d$ is non-negative integer integer it can be used to ﬁnd $ω$ using $\begin{aligned} ω & = \sum_{c \in Γ} (s (c) [\sqrt{r}]_{c} + \frac{α_{c}^{s (c)}}{x - c}) + s (\infty) [\sqrt{r}]_{\infty} \end{aligned}$

The above reduces to $\begin{aligned} ω & = (-) [\sqrt{r}]_{\infty} \\ = 0 + (-) (\frac{i \sqrt{3}}{2}) \\ = - \frac{i \sqrt{3}}{2} \end{aligned}$

Now that $ω$ is determined, the next step is ﬁnd a corresponding minimal polynomial $p (x)$ of degree $d = 0$ to solve the ode. The polynomial $p (x)$ needs to satisfy the equation $\begin{aligned} (1) & p^{″} + 2 ω p^{'} + (ω^{'} + ω^{2} - r) p & = 0 \end{aligned}$

Since $d = 0$ then let $p (x) = 1$ . Substituting this in the above gives $\begin{aligned} (0) + 2 (- \frac{i \sqrt{3}}{2}) (0) + ((0) + {(- \frac{i \sqrt{3}}{2})}^{2} - (- \frac{3}{4})) & = 0 \\ 0 & = 0 \end{aligned}$

The equation is satisﬁed. Therefore the ﬁrst solution to the ode $z^{″} = r z$ is $\begin{aligned} z (x) & = p e^{\int ω d x} \\ = e^{\int - \frac{i \sqrt{3}}{2} d x} \\ = e^{- \frac{i \sqrt{3} x}{2}} \end{aligned}$

The ﬁrst solution to the original ode in $y$ is now found from (using $A = 1, B = 1$ ) $\begin{aligned} y_{1} & = z e^{\int - \frac{1}{2} \frac{B}{A} d x} \\ = z e^{- \int \frac{1}{2} d x} \\ = z e^{- \frac{x}{2}} \\ = e^{- \frac{i \sqrt{3} x}{2}} (e^{- \frac{x}{2}}) \\ = e^{- \frac{x (1 + I \sqrt{3})}{2}} \end{aligned}$

The second solution $y_{2}$ to the original ode is found using reduction of order $\begin{aligned} y_{2} & = y_{1} \int \frac{e^{\int - \frac{B}{A} d x}}{y_{1}^{2}} d x \\ = y_{1} \int \frac{e^{- \int d x}}{{(y_{1})}^{2}} d x \\ = y_{1} \int \frac{e^{- x}}{{(y_{1})}^{2}} d x \\ = (e^{- \frac{x (1 + I \sqrt{3})}{2}}) (- \frac{i \sqrt{3} e^{i \sqrt{3} x}}{3}) \\ = - \frac{I}{3} e^{\frac{x (I \sqrt{3} - 1)}{2}} \sqrt{3} \end{aligned}$

Therefore the general solution is $\begin{aligned} y & = c_{1} y_{1} + c_{2} y_{2} \\ = c_{1} (e^{- \frac{x (1 + i \sqrt{3})}{2}}) + c_{2} (- \frac{i e^{\frac{x (i \sqrt{3} - 1)}{2}} \sqrt{3}}{3}) \end{aligned}$

Using Euler’s formula the above can be simpliﬁed to the standard looking solution $\begin{aligned} y (x) & = e^{- \frac{x}{2}} (C_{1} \sin (\frac{\sqrt{3} x}{2}) + C_{2} \cos (\frac{\sqrt{3} x}{2})) \end{aligned}$

3.2 case two

3.2.1 Example 1

Given the ode $\begin{aligned} 2 x^{2} y^{″} - x y^{'} + (1 + x) y & = 0 \end{aligned}$

Converting it $y^{″} + a y^{'} + b y = 0$ gives $\begin{aligned} y^{″} - \frac{1}{2 x} y^{'} + \frac{1 + x}{2 x^{2}} y & = 0 \end{aligned}$

Where $a = - \frac{1}{2 x}, b = \frac{1 + x}{2 x^{2}}$ . Applying the transformation $z = y e^{\frac{1}{2} \int a d x}$ gives $z^{″} = r z$ and $r = \frac{1}{4} a^{2} + \frac{1}{2} a^{'} - b$ . Therefore $\begin{aligned} r & = \frac{s}{t} \\ = \frac{- 3 - 8 x}{16 x^{2}} \end{aligned}$

There is one pole at $x = 0$ of order $2$ and $O (\infty) = \deg (t) - \deg (s) = 2 - 1 = 1$ . Table 1 shows that necessary conditions for only case two are satisﬁed. Since pole is order $2$ , the set $E_{0} = {2, 2 + 2 \sqrt{1 + 4 b}, 2 - 2 \sqrt{1 + 4 b}}$ where $b$ is the coeﬃcient of $\frac{1}{(x)^{2}}$ in the partial fraction decomposition of $r$ given by $\begin{aligned} r & = - \frac{3}{16 x^{2}} - \frac{1}{2 x} \end{aligned}$

Therefore $b = - \frac{3}{16}$ . Hence $\begin{aligned} E_{0} & = {2, 2 + 2 \sqrt{1 - 4 \frac{3}{16}}, 2 - 2 \sqrt{1 - 4 \frac{3}{16}}} \\ = {2, 3, 1} \end{aligned}$

Since $O (\infty) = 1 < 2$ then $E_{\infty} = O (\infty) = {1}$ . This completes step $1$ of case two. Step $2$ is used to determine a non-negative integer $d$ . Using $e_{\infty} = 1$ gives $\begin{aligned} d & = \frac{1}{2} (e_{\infty} - \sum_{c \in Γ} e_{c}) \end{aligned}$

There is only one pole, so the sum contains only one term. There are 3 possible combinations to try, using either $e_{0} = 2$ , $e_{0} = 3$ or $e_{0} = 1$ . Therefore $\begin{aligned} d & = \frac{1}{2} (e_{\infty} - (e_{0})) & = \frac{1}{2} (1 - 2) & = - \frac{1}{2} \\ = \frac{1}{2} (e_{\infty} - (e_{0})) & = \frac{1}{2} (1 - 3) & = - 1 \\ = \frac{1}{2} (e_{\infty} - (e_{0})) & = \frac{1}{2} (1 - 1) & = 0 \end{aligned}$

The above shows that only the family ${e_{\infty} = 1, e_{0} = 1}$ generated non-negative $d = 0$ . $θ$ is now found. In the following sum, only $e_{c}$ retained from the above are used. In this example, this will be $e_{0} = 1$ since it is the member of $E_{0}$ which generated non-negative integer. If there were more than one $e_{i}$ found, then each would be tried at time. $\begin{aligned} θ & = \frac{1}{2} \sum_{c \in Γ} \frac{e_{c}}{x - c} \\ = \frac{1}{2} (\frac{e_{0}}{x - 0}) \\ = \frac{1}{2 x} \end{aligned}$

This completes step $2$ . Step $3$ ﬁnds polynomial $p (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + x^{d}$ of degree $d$ . Since $d = 0$ then $p (x) = 1$ . This polynomial has to satisfy the following $\begin{array}{r} p^{‴} + 3 θ p^{″} (3 θ^{2} + 3 θ^{'} - 4 r) p^{'} + (θ^{″} + 3 θ θ^{'} + θ^{3} - 4 r θ - 2 r^{'}) p = 0 \end{array}$

Substituting $p = 1, θ = \frac{1}{2 x}$ into the above and simplifying gives $\begin{aligned} 0 & = 0 \end{aligned}$

Since $p (x) = 1$ is veriﬁed, then $\begin{array}{r} ϕ = θ + \frac{p^{'}}{p} \\ = \frac{1}{2 x} \end{array}$

Next, $ω$ solution is found using $\begin{aligned} ω^{2} - ϕ ω + (\frac{1}{2} ϕ^{'} + \frac{1}{2} ϕ^{2} - r) & = 0 \end{aligned}$

Substituting the values for $ϕ$ and $r$ into the above gives $\begin{array}{r} w^{2} - \frac{w}{2 x} + \frac{1 + 8 x}{16 x^{2}} = 0 \end{array}$

Solving for $ω$ gives two roots, either one can be used. Using $\begin{aligned} ω & = \frac{1 + 2 \sqrt{2} \sqrt{- x}}{4 x} \end{aligned}$

Therefore the ﬁrst solution to the ode $z^{″} = r z$ is $\begin{aligned} z (x) & = e^{\int ω d x} \\ = e^{\int \frac{1 + 2 \sqrt{2} \sqrt{- x}}{4 x} d x} \\ = x^{\frac{1}{4}} e^{\sqrt{2} \sqrt{- x}} \end{aligned}$

The ﬁrst solution to the original ode in $y$ is now found from $\begin{aligned} y_{1} & = z e^{\int - \frac{1}{2} a d x} \\ = z e^{- \int \frac{1}{2} \frac{- x}{2 x^{2}} d x} \\ = z e^{\frac{\ln (x)}{4}} \\ = (x^{\frac{1}{4}} e^{\sqrt{2} \sqrt{- x}}) x^{\frac{1}{4}} \\ = \sqrt{x} e^{\sqrt{2} \sqrt{- x}} \end{aligned}$

The second solution $y_{2}$ to the original ode can be found using reduction of order.

3.2.2 Example 2

This is an ode in which the necessary conditions for all three cases are satisﬁed, but solved using case two to illustrate the algorithm. $\begin{array}{r} (1 - x) x^{2} y^{″} + (5 x - 4) x y^{'} + (6 - 9 x) y = 0 \end{array}$

Converting it $y^{″} + a y^{'} + b y = 0$ gives $\begin{array}{r} y^{″} + \frac{5 x - 4}{(1 - x) x} y^{'} + \frac{6 - 9 x}{(1 - x) x^{2}} y = 0 \end{array}$

Where $a = \frac{5 x - 4}{(1 - x) x}, b = \frac{6 - 9 x}{(1 - x) x^{2}}$ . Applying the transformation $z = y e^{\frac{1}{2} \int a d x}$ gives $z^{″} = r z$ where $r = \frac{1}{4} a^{2} + \frac{1}{2} a^{'} - b$ . Therefore $\begin{aligned} r & = \frac{s}{t} \\ = \frac{4 - x}{4 x (x - 1)^{2}} \end{aligned}$

There is a pole at $x = 0$ of order $1$ and a pole at $x = 1$ of order $2$ . Since there is no odd order pole larger than $2$ and the order at $\infty$ is $2$ then the necessary conditions for case one are satisﬁed. Since there is a pole of order $2$ then the necessary conditions for case two are also satisﬁed. Since pole order is not larger than $2$ and the order at $\infty$ is $2$ then the necessary conditions for case three are also satisﬁed. Any one of the three cases algorithm could be used to solve this, but here case two will be used for illustration.

The pole of order $1$ at $x = 0$ gives $E_{0} = {4}$ and the pole of order $2$ at $x = 1$ gives $E_{1} = {2, 2 + 2 \sqrt{1 + 4 b}, 2 - 2 \sqrt{1 + 4 b}}$ where $b$ is the coeﬃcient of $\frac{1}{(x - 1)^{2}}$ in the partial fraction decomposition of $r$ . The partial fractions decomposition of $r$ is $\begin{array}{r} r = \frac{3}{4 {(x - 1)}^{2}} - \frac{1}{x - 1} + \frac{1}{x} \end{array}$

The above shows that $b = \frac{3}{4}$ , therefore $E_{1} = {- 2, 2, 6}$ .

$O (\infty) = 2$ therefore $E_{\infty} = {2, 2 + 2 \sqrt{1 + 4 b}, 2 - 2 \sqrt{1 + 4 b}}$ where $b = \frac{lcoef (s)}{lcoeff (t)}$ where $r = \frac{s}{t}$ . This gives $b = - \frac{1}{4}$ , hence $\begin{aligned} E_{\infty} & = {2, 2 + 2 \sqrt{1 + 4 b}, 2 - 2 \sqrt{1 + 4 b}} \\ = {2, 2 + 2 \sqrt{1 - 4 \frac{1}{4}}, 2 - 2 \sqrt{1 - 4 - \frac{1}{4}}} \\ = {2, 2, 2} \\ = {2} \end{aligned}$

The following table summarizes step $1$ results.


pole $c$ location	pole order	$E_{c}$

$0$	$1$	${4}$

$1$	$2$	${- 2, 2, 6}$

Table 4: First step, case two.

E_{c}

set information


Order of $r$ at $\infty$	$E_{\infty}$

$2$	${2}$

Table 5: First step, case two.

O (\infty)

information

The above completes step 1 for case two. Step 2 searches for a non-negative integer $d$ using $\begin{aligned} d & = \frac{1}{2} (e_{\infty} - \sum_{c \in Γ} e_{c}) \end{aligned}$

Where in the above $e_{c} \in E_{c}$ , $e_{\infty} \in E_{\infty}$ were found in the ﬁrst step. The following are the possible combinations to use $\begin{aligned} d & = \frac{1}{2} (2 - (4 - 2)) & = 0 \\ = \frac{1}{2} (2 - (4 + 2)) & = - 2 \\ = \frac{1}{2} (2 - (4 + 6)) & = - 4 \end{aligned}$

The above shows that $E_{c} = {e_{0} = 4, e_{1} = - 2}$ are the family of values to use and all other values are discarded.

The following rational function $θ$ is determined using $\begin{aligned} θ & = \frac{1}{2} \sum_{c \in Γ} \frac{e_{c}}{x - c} \\ = \frac{1}{2} (\frac{4}{(x - (0))} + \frac{- 2}{(x - (1))}) \\ = \frac{2}{x} - \frac{1}{x - 1} \end{aligned}$

The algorithm now searches for a monic polynomial $p (x)$ of degree $d = 0$ such that $\begin{aligned} p^{‴} + 3 θ p^{″} + (3 θ^{2} + 3 θ^{'} - 4 r) p^{'} + (θ^{″} + 3 θ θ^{'} + θ^{3} - 4 r θ - 2 r^{'}) p & = 0 \end{aligned}$

Since $d = 0$ , then $p (x) = 1$ . Substituting the values found in step $2$ in the above equation and simplifying gives $\begin{array}{r} 0 = 0 \end{array}$

Hence $p (x) = 1$ can be used. Let $\begin{aligned} ϕ & = θ + \frac{p^{'}}{p} \\ = \frac{2}{x} - \frac{1}{x - 1} \end{aligned}$

And $ω$ be the solution of $\begin{aligned} ω^{2} - ϕ ω + (\frac{1}{2} ϕ^{'} + \frac{1}{2} ϕ^{2} - r) & = 0 \end{aligned}$

Substituting the values for $ϕ$ and $r$ into the above equation gives $\begin{aligned} w^{2} - (\frac{2}{x} - \frac{1}{x - 1}) w + \frac{{(x - 2)}^{2}}{4 {(x - 1)}^{2} x^{2}} & = 0 \end{aligned}$

Solving for $ω$ gives two roots, either one can be used. Using $\begin{aligned} ω & = \frac{x - 2}{2 (x - 1) x} \end{aligned}$

Therefore to $z^{″} = r z$ is $\begin{aligned} z & = e^{\int ω d x} \\ = e^{\int \frac{x - 2}{2 (x - 1) x} d x} \\ = \frac{x}{\sqrt{x - 1}} \end{aligned}$

The ﬁrst solution to the original ode in $y$ is now found from $\begin{aligned} y_{1} & = z e^{- \int \frac{1}{2} a d x} \\ = z e^{- \int \frac{1}{2} \frac{5 x^{2} - 4 x}{- x^{3} + x^{2}} d x} \\ = z e^{2 \ln (x) + \frac{\ln (x - 1)}{2}} \\ = z (x^{2} \sqrt{x - 1}) \\ = \frac{x}{\sqrt{x - 1}} (x^{2} \sqrt{x - 1}) \\ = x^{3} \end{aligned}$

The second solution $y_{2}$ to the original ode is found using reduction of order.

3.3 case three

3.3.1 Example 1

This is the same ode used in second example above for case two as the necessary conditions for case three are also satisﬁed. $\begin{array}{r} (1 - x) x^{2} y^{″} + (5 x - 4) x y^{'} + (6 - 9 x) y = 0 \end{array}$

As shown earlier, this ode is transformed to $z^{″} = r z$ where $\begin{aligned} r & = \frac{s}{t} \\ = \frac{4 - x}{4 x (x - 1)^{2}} \end{aligned}$

There is a pole at $x = 0$ of order $1$ and a pole at $x = 1$ of order $2$ . For the pole of order $1$ at $x = 0$ , $E_{0} = {12}$ . For the pole of order $2$ at $x = 1$ $\begin{aligned} (1) & E_{1} & = {6 + \frac{12 k}{n} \sqrt{1 + 4 b}} for k & = - \frac{n}{2} \dots \frac{n}{2} \end{aligned}$

Where $k$ is incremented by $1$ each time, and $n$ is any of ${4, 6, 12}$ and $b$ is the coeﬃcient of $\frac{1}{(x - 1)^{2}}$ in the partial fraction decomposition of $r$ which is $\begin{array}{r} r = \frac{3}{4 {(x - 1)}^{2}} - \frac{1}{x - 1} + \frac{1}{x} \end{array}$

The above shows that $b = \frac{3}{4}$ . Starting with $n = 4$ (if this $n$ produces no solution then $n = 6, 12$ will be tried also). Equation (1) now becomes $\begin{aligned} E_{1} & = {6 + \frac{12 k}{4} \sqrt{1 + 4 (\frac{3}{4})}} for k & = - 2 \dots 2 \end{aligned}$

Which simpliﬁes to $\begin{aligned} E_{1} & = {6 + 6 k} for k & = - 2 \dots 2 \\ = {- 6, 0, 6, 12, 18} \end{aligned}$

$E_{\infty}$ is found using (1) but with diﬀerent $b$ . In this case $b$ is given by $b = \frac{lcoef (s)}{lcoeff (t)}$ where $r = \frac{s}{t}$ . $lcoef (s)$ is the leading coeﬃcient of $s$ and $lcoef (t)$ is the leading coeﬃcient of $t$ . Since $r = \frac{- x + 4}{4 x^{3} - 8 x^{2} + 4 x}$ then $b = - \frac{1}{4}$ . Equation (1) becomes $\begin{aligned} E_{\infty} & = {6 + \frac{12 k}{4} \sqrt{1 - 4 (\frac{1}{4})}} for k & = - 2 \dots 2 \end{aligned}$

This simpliﬁes to $\begin{aligned} E_{\infty} & = {6} \end{aligned}$

The following table summarizes step $1$ results using $n = 4$ .


pole $c$ location	pole order	$E_{c}$

$0$	$1$	${12}$

$1$	$2$	${- 6, 0, 6, 12, 18}$

Table 6: First step, case three using

n = 4

E_{c}

set information


Order of $r$ at $\infty$	$E_{\infty}$

$2$	${6}$

Table 7: First step, case three using

n = 4

O (\infty)

information

The next step is to determine a non negative integer $d$ using $\begin{aligned} d & = \frac{n}{12} (e_{\infty} - \sum_{c \in Γ} e_{c}) \end{aligned}$

Where in the above $e_{c}$ is a distinct element from each corresponding $E_{c}$ . This means all possible tuples ${e_{c_{1}}, e_{c_{2}}, \dots, e_{c_{n}}}$ are tried in the sum above, where $e_{c_{i}}$ is one element of each $E_{c}$ found earlier.

This results in the following values for $d$ using $n = 4$ and $e_{\infty} = 6$ . $\begin{aligned} d & = \frac{1}{3} (6 - (12 - 6)) & = 0 \\ = \frac{1}{3} (6 - (12 + 0)) & = - 2 \\ = \frac{1}{3} (6 - (12 + 6)) & = - 4 \\ = \frac{1}{3} (6 - (12 + 12)) & = - 6 \\ = \frac{1}{3} (6 - (12 + 18)) & = - 8 \end{aligned}$

Therefore only the ﬁrst case using $e_{0} = 12, e_{1} = - 6$ generated non-negative integer $d$ . The following rational function is now formed $\begin{aligned} θ & = \frac{n}{12} \sum_{c \in Γ} \frac{e_{c}}{x - c} \\ = \frac{4}{12} (\frac{12}{(x - (0))} + \frac{- 6}{(x - (1))}) \\ = \frac{2 x - 4}{(x - 1) x} \end{aligned}$

And $\begin{aligned} S & = \prod_{c \in Γ} (x - c) \\ = (x - 0) (x - 1) \\ = x (x - 1) \end{aligned}$

This completes the step 2 of the algorithm.

Since the degree $d = 0$ , then $p (x) = 1$ . Now $P_{i} (x)$ polynomials are generated using $\begin{aligned} P_{n} & = - p (x) \\ P_{i - 1} & = - S p_{i}^{'} + ((n - i) S^{'} - S θ) P_{i} - (n - 1) (i + 1) S^{2} r P_{i + 1} i = n, n - 1, \dots, 0 \end{aligned}$

These generate the following set $\begin{aligned} P_{4} & = - 1 \\ P_{3} & = 2 x - 4 \\ P_{2} & = - 3 {(x - 2)}^{2} \\ P_{1} & = 3 {(x - 2)}^{3} \\ P_{0} & = - \frac{3 {(x - 2)}^{4}}{2} \\ P_{- 1} & = 0 \end{aligned}$

There is nothing to solve for from the last equation $P_{- 1} = 0$ as $p (x) = 1$ is already known because the degree $d$ was zero and hence there are no coeﬃcients $a_{i}$ to solve for.

Next $ω$ is determined as the solution to the following equation using $n = 4$ . $\begin{aligned} \sum_{i = 0}^{n} S^{i} \frac{P_{i}}{(n - i)!} ω^{i} & = 0 \\ \frac{P_{0}}{4!} + \frac{S P_{1}}{3!} ω + \frac{S^{2} P_{2}}{2!} ω^{2} + \frac{S^{3} P_{3}}{1!} ω^{3} + + \frac{S^{4} P_{4}}{0!} ω^{4} & = 0 \\ - \frac{1}{16} {(2 ω x^{2} - 2 x ω - x + 2)}^{4} & = 0 \end{aligned}$

Solving the above and using any one of the roots results in $\begin{array}{r} ω = \frac{1}{2 x (x - 1)} (x - 2) \end{array}$

The above $ω$ is used to ﬁnd a solution to $z^{″} = r z$ from $\begin{aligned} z & = e^{\int ω d x} \\ = e^{\int \frac{x - 2}{2 x (x - 1)} d x} \\ = \frac{x}{\sqrt{x - 1}} \end{aligned}$

Therefore one solution to the original ode in $y$ is $\begin{aligned} y & = z e^{\int - \frac{1}{2} a d x} \\ = z e^{- \int \frac{1}{2} \frac{5 x^{2} - 4 x}{- x^{3} + x^{2}} d x} \\ = z e^{2 \ln (x) + \frac{\ln (x - 1)}{2}} \\ = (\frac{x}{\sqrt{x - 1}}) (x^{2} \sqrt{x - 1}) \end{aligned}$

The second solution to the original ode is found using reduction of order. This completes the solution using case $3$ for degree $n = 4$ of $ω$ .

3.3.2 Example 2

The ode is $\begin{aligned} x^{2} (1 + x) y^{″} + x (2 x + 1) y^{″} - (4 + 6 x) y & = 0 \end{aligned}$

Converting it $y^{″} + a y^{'} + b y = 0$ gives $\begin{array}{r} y^{″} + \frac{2 x + 1}{x (1 + x)} y^{'} - \frac{4 + 6 x}{x^{2} (1 + x)} y = 0 \end{array}$

Where $a = \frac{2 x + 1}{x (1 + x)}, b = - \frac{4 + 6 x}{x^{2} (1 + x)}$ . Applying the transformation $z = y e^{\frac{1}{2} \int a d x}$ results in $z^{″} = r z$ where $r = \frac{1}{4} a^{2} + \frac{1}{2} a^{'} - b$ where $\begin{aligned} r & = \frac{s}{t} \\ = \frac{24 x^{2} + 40 x + 15}{4 {(x (x + 1))}^{2}} \end{aligned}$

There is a pole at $x = 0$ of order $2$ and a pole at $x = - 1$ of order $2$ . Since there is no odd order pole larger than $2$ and the order at $\infty$ is $2$ then the necessary conditions for case one are satisﬁed. Since there is a pole of order $2$ then necessary conditions for case two are also satisﬁed. Since pole order is not larger than $2$ and the order at $\infty$ is $2$ then the necessary conditions for case three are also satisﬁed. This is now solved as case three for illustration.

Starting with $n = 4$ , and for the pole of order $2$ at $x = - 1$ $\begin{aligned} E_{- 1} & = {6 + \frac{12 k}{n} \sqrt{1 + 4 b}} for k & = - \frac{n}{2} \dots \frac{n}{2} \end{aligned}$

This simpliﬁes to $\begin{aligned} (1) & E_{- 1} & = {6 + 3 k \sqrt{1 + 4 b}} for k & = - 2 \dots 2 \end{aligned}$

$b$ is the coeﬃcient of $\frac{1}{{(1 + x)}^{2}}$ in the partial fractions decomposition of $r$ given by $\begin{array}{r} r = - \frac{1}{4 {(1 + x)}^{2}} + \frac{15}{4 x^{2}} - \frac{5}{2 (1 + x)} + \frac{5}{2 x} \end{array}$

The above shows that $b = - \frac{1}{4}$ . Equation (1) becomes $\begin{aligned} E_{- 1} & = {6 + 3 k \sqrt{1 - 4 (\frac{1}{4})}} for k & = - 2 \dots 2 \\ = {6} \end{aligned}$

For the pole at $x = 0$ of order $2$ , $b$ is the coeﬃcient of $\frac{1}{x^{2}}$ in the above partial fractions decomposition of $r$ . This shows that $b = \frac{15}{4}$ . Hence $\begin{aligned} E_{0} & = {6 + 3 k \sqrt{1 + 4 (\frac{15}{4})}} for k & = - 2 \dots 2 \\ = {- 18, - 6, 6, 18, 30} \end{aligned}$

$E_{\infty}$ is found using equation (1) but with diﬀerent $b$ . In this case $b$ is given by $b = \frac{lcoef (s)}{lcoeff (t)}$ where $r = \frac{s}{t}$ . $lcoef (s)$ is the leading coeﬃcient of $s$ and $lcoef (t)$ is the leading coeﬃcient of $t$ . Since $r = \frac{24 x^{2} + 40 x + 15}{4 x^{4} + 8 x^{3} + 4 x^{2}}$ then $b = 6$ . Equation (1) becomes $\begin{aligned} E_{\infty} & = {6 + 3 k \sqrt{1 + 4 (6)}} for k & = - 2 \dots 2 \\ = {- 24, - 9, 6, 21, 36} \end{aligned}$

The following table summarizes step $1$ results using $n = 4$ .


pole $c$ location	pole order	$E_{c}$

$- 1$	$2$	${6}$

$0$	$2$	${- 18, - 6, 6, 18, 30}$

Table 8: First step, case three using

n = 4

E_{c}

set information


Order of $r$ at $\infty$	$E_{\infty}$

$2$	${- 24, - 9, 6, 21, 36}$

Table 9: First step, case three using

n = 4

O (\infty)

information

The next step is to determine a non negative integer $d$ using $\begin{aligned} d & = \frac{n}{12} (e_{\infty} - \sum_{c \in Γ} e_{c}) \end{aligned}$

This results in the following values for $d$ using $n = 4$ . $\begin{aligned} d & = \frac{1}{3} (- 24 - (6 - 18)) & = - 4 \\ = \frac{1}{3} (- 24 - (6 + 6)) & = - 12 \\ = \frac{1}{3} (- 24 - (6 - 6)) & = - 8 \\ = \frac{1}{3} (- 24 - (6 + 18)) & = - 16 \\ = \frac{1}{3} (- 24 - (6 + 30)) & = - 20 \\ = \frac{1}{3} (- 9 - (6 - 18)) & = 1 \\ = \frac{1}{3} (- 9 - (6 + 6)) & = - 7 \\ = \frac{1}{3} (- 9 - (6 - 6)) & = - 3 \\ = \frac{1}{3} (- 9 - (6 + 18)) & = - 11 \\ = \frac{1}{3} (- 9 - (6 + 30)) & = - 15 \\ = \frac{1}{3} (6 - (6 - 18)) & = 6 \\ = \frac{1}{3} (6 - (6 + 6)) & = - 2 \\ = \frac{1}{3} (6 - (6 - 6)) & = 2 \\ = \frac{1}{3} (6 - (6 + 18)) & = - 6 \\ = \frac{1}{3} (6 - (6 + 30)) & = - 10 \\ = \frac{1}{3} (21 - (6 - 18)) & = 11 \\ = \frac{1}{3} (21 - (6 + 6)) & = 3 \\ = \frac{1}{3} (21 - (6 - 6)) & = 7 \\ = \frac{1}{3} (21 - (6 + 18)) & = - 1 \\ = \frac{1}{3} (21 - (6 + 30)) & = - 5 \\ = \frac{1}{3} (36 - (6 - 18)) & = 16 \\ = \frac{1}{3} (36 - (6 + 6)) & = 8 \\ = \frac{1}{3} (36 - (6 - 6)) & = 12 \\ = \frac{1}{3} (36 - (6 + 18)) & = 4 \\ = \frac{1}{3} (36 - (6 + 30)) & = 0 \end{aligned}$

From the above, the following families all produce non-negative $d$

$e_{\infty} = - 9$	$e_{- 1} = 6$	$e_{0} = - 18$	$d = 1$
$e_{\infty} = 6$	$e_{- 1} = 6$	$e_{0} = - 18$	$d = 6$
$e_{\infty} = 6$	$e_{- 1} = 6$	$e_{0} = - 6$	$d = 2$
$e_{\infty} = 21$	$e_{- 1} = 6$	$e_{0} = - 18$	$d = 11$
$e_{\infty} = 21$	$e_{- 1} = 6$	$e_{0} = 6$	$d = 3$
$e_{\infty} = 21$	$e_{- 1} = 6$	$e_{0} = - 6$	$d = 7$
$e_{\infty} = 36$	$e_{- 1} = 6$	$e_{0} = - 18$	$d = 16$
$e_{\infty} = 36$	$e_{- 1} = 6$	$e_{0} = 6$	$d = 8$
$e_{\infty} = 36$	$e_{- 1} = 6$	$e_{0} = - 6$	$d = 12$
$e_{\infty} = 36$	$e_{- 1} = 6$	$e_{0} = 18$	$d = 4$
$e_{\infty} = 36$	$e_{- 1} = 6$	$e_{0} = 30$	$d = 0$

Starting with the smallest $d = 0$ as that is the least computationally expensive with its corresponding $e_{- 1} = 6, e_{0} = 30, e_{\infty} = 36$ , the following rational function is now formed $\begin{aligned} θ & = \frac{n}{12} \sum_{c \in Γ} \frac{e_{c}}{x - c} \\ = \frac{4}{12} (\frac{6}{x + 1} + \frac{30}{x}) \\ = \frac{12 x + 10}{x (1 + x)} \end{aligned}$

And $\begin{aligned} S & = \prod_{c \in Γ} (x - c) \\ = (x + 1) x \end{aligned}$

This completes the step 2 of the algorithm.

Polynomial $p (x)$ is now determined. Since the degree of the polynomial is $d = 0$ , then $\begin{array}{r} p (x) = 1 \end{array}$

The $P_{i} (x)$ polynomials are generated using $\begin{aligned} P_{n} & = - p (x) \\ P_{i - 1} & = - S p_{i}^{'} + ((n - i) S^{'} - S θ) P_{i} - (n - 1) (i + 1) S^{2} r P_{i + 1} i = n, n - 1, \dots, 0 \end{aligned}$

The above results in the following set $\begin{aligned} P_{4} & = - 1 \\ P_{3} & = 12 x + 10 \\ P_{2} & = - 3 {(6 x + 5)}^{2} \\ P_{1} & = 3 {(6 x + 5)}^{3} \\ P_{0} & = - \frac{3 {(6 x + 5)}^{4}}{2} \\ P_{- 1} & = 0 \end{aligned}$

There is coeﬃcient for $p (x)$ to solve for from the last equation $P_{- 1} = 0$ as $p (x) = 1$ is already known because the degree $d$ is zero. $ω$ is now determined as the solution to the following equation, using $n = 4$ $\begin{aligned} \sum_{i = 0}^{n} S^{i} \frac{P_{i}}{(n - i)!} ω^{i} & = 0 \\ \frac{P_{0}}{4!} + \frac{S P_{1}}{3!} ω + \frac{S^{2} P_{2}}{2!} ω^{2} + \frac{S^{3} P_{3}}{1!} ω^{3} + + \frac{S^{4} P_{4}}{0!} ω^{4} & = 0 \\ - \frac{1}{16} {(2 ω x^{2} + 2 x ω - 6 x - 5)}^{4} & = 0 \end{aligned}$

Solving the above and using any one of the roots gives $\begin{array}{r} ω = \frac{x - 2}{2 x (x - 1)} \end{array}$

This $ω$ is used to ﬁnd a solution to $z^{″} = r z$ from $\begin{aligned} z & = e^{\int ω d x} \\ = e^{\int \frac{6 x + 5}{2 x (1 + x)} d x} \\ = x^{\frac{5}{2}} \sqrt{1 + x} \end{aligned}$

The ﬁrst solution to the original ode in $y$ is found from $\begin{aligned} y & = z e^{\int - \frac{1}{2} a d x} \\ = z e^{- \int \frac{1}{2} \frac{2 x^{2} + x}{x^{2} (1 + x)} d x} \\ = z e^{- \frac{\ln (x (1 + x))}{2}} \\ = (x^{\frac{5}{2}} \sqrt{1 + x}) (\frac{1}{\sqrt{x (1 + x)}}) \end{aligned}$

Which simpliﬁes to $\begin{aligned} y & = x^{2} \end{aligned}$

The second solution to the original ode is found using reduction of order.

4 Conclusion

Detailed description of the Kovacic algorithm with worked out examples were given. All three cases of the Kovacic algorithm were implemented using object oriented design in Maple. The software was then used to analyze over $3000$ diﬀerential equations. The results showed that case one and two combined provided coverage for $99.9$ % of the ode’s with $97.36$ % of the ode’s solved using case one algorithm and $2.54$ % solved using case $2$ algorithm with only $0.1$ % requiring case $3$ . Not a single ode was found that required the use of case three with $n = 6$ or $n = 12$ .

One restriction found on the use of the algorithm is that it requires an ode with its coeﬃcients being numerical and not symbolic. This is because the algorithm has to decide in step $2$ if $d$ (the degree of polynomial $p (x)$ ) is non-negative integer or not in order to continue to step $3$ . If some of the ode coeﬃcients were symbolic, it will not be able to decide on this (without additional assumptions provided). Therefore this algorithm works best with ode’s having its coeﬃcients given with numerical values.

5 Appendix

5.1 Instructions and examples using the Kovacic package

The Kovacic class is included in the ﬁle KOV.mpl and the Kovacic testsuite module is in the ﬁle kovacic_tester.mpl. These two ﬁles accompany the arXiv version of this paper.

To use these, download these two ﬁles to some directory at your computer. For example, on windows, assuming the ﬁles were downloaded to c:/my_folder/, then now start Maple and type

read "c:/my_folder/KOV.mpl" 
read "c:/my_folder/kovacic_tester.mpl"

The above will load the kovacic_class and the testsuite module. Once the above is successfully completed, then to solve an ode the command is

ode  := diff(y(x),x$2)+diff(y(x),x)+y(x)=0; 
o    := Object(kovacic_class,ode,y(x)); #create the object 
sol  := o:-dsolve();

The above command will automatically try all the cases that have been detected one by one until a solution is found. If no solution is found, it returns FAIL. To verify the solution, the command is

if sol<>FAIL then 
   odetest(sol,ode); 
fi;

Which returns $0$ if the solution is correct.

A note on the type of ode’s supported: it is recommenced to use only ode’s with numeric coeﬃcients and not symbolic coeﬃcients. This is because the Kovacic algorithm needs to decide if the degree $d$ of the polynomial $p (x)$ is non-negative or not. If some of the coeﬃcients are purely symbolic, then it can fail to decide this. An example of this is given in the original Kovacic paper as example 2 on page 14, which is to solve the Bessel ode $y^{″} = (\frac{4 n^{2} - 1}{4 x^{2}} - 1) y$ . This will now return FAIL since the algorithm can not decide if $d$ is non-negative integer without knowing any assumptions or having numerical value for $n$ . Replacing $n$ by any half odd integer, then it can solve it as follows

ode :=diff(diff(y(x),x),x)=((4*n^2-1)/(4*x^2)-1)*y(x); 
n   :=-3/2; 
o   := Object(kovacic_class,ode,y(x)); 
sol := o:-dsolve();

To solve an ode using speciﬁc case number, say case $2$ , the command is

ode   := ...; 
o     := Object(kovacic_class,ode,y(x)); 
sol   := o:-dsolve_case(2);

If the ode happened to satisfy cases $1$ and $2$ for an example, then the above command will only use case $2$ to solve it and will skip case $1$ . If the command o:-dsolve() was used instead, then the ode will be solved using case $1$ instead as that is the ﬁrst one tried. Case $2$ will only be tried is no solution is found using case $1$ .

The object created above, named as “o”, has additional public methods that can be invoked. The following is description of all public methods available.

o:-get_y_ode() This returns the original ode.
o:-get_z_ode() This returns the ode solved by Kovacic algorithm which is $z^{″} = r z$ .
o:-get_r() This returns $r$ only.
o:-get_poles() This returns list of the poles of $r$ . It has the format

[ [pole location,pole order],[pole location,pole order], ...]

If there are no poles, then the empty list [] is returned.
o:-get_order_at_infinity() This returns the order of $r$ at inﬁnity.
o:-get_possible_cases() This returns list of possible Kovacic cases detected which can be [1], [2], [1,2], [1,2,3]. If no Kovacic cases are found, then the empty list [] is returned.
o:-get_case_used() This returns the actual case number used if solution to the ode was successful. This can be $1, 2$ or $3$ . If no solution is found after trying all cases whose conditions were satisﬁed, then $- 1$ is returned.
o:-get_n_case_3() This return $n$ used when case $3$ was used to solve the ode. This can be $4, 6$ or $12$ . If case $3$ was not used, or no solution is found, then $- 1$ is returned.

To run the full testsuite of $3000$ ode’s that comes with the package, the command is

kovacic_tester:-unit_test_main_api(); 
                   "Test ", 6735, " PASSED " 
                   "Test ", 6736, " PASSED " 
                   "Test ", 6737, " PASSED " 
                             . 
                             . 
                   "Test ", 7579, " PASSED " 
                   "Test ", 7580, " PASSED " 
                   "Test ", 7581, " PASSED "

To run testsuite using speciﬁc cases only, the commands are

kovacic_tester:-unit_test_case_1(); 
kovacic_tester:-unit_test_case_2(); 
kovacic_tester:-unit_test_case_3();

5.2 Source code

#-------------------------------------------------------------------- 
#FILE NAME  :   KOV.mpl 
# 
# Copyright (C) 2022, Nasser M. Abbasi. All rights reserved 
# email: nma@12000.org 
# 
# Free software to use and modify in anyway as long as the above 
# copyright notice remains attached in the file 
# 
# Change history 
#----------------- 
# Oct 27, 2022.     Initial version 
# 
# Note that the latest version and any updates can alaways be obtained 
# from the author web site at 
# http://12000.org/my_notes/my_paper_on_kovacic/paper.htm 
# 
# Any problems found in the software please report so I can correct. 
# This implementation was done using Maple 2022.2 on windows 10. 
# 
#-------------------------------------------------------------------- 
 
#-------------------------------------------------------------------- 
# An Object oriented implementation of the Kovacic algortithm 
# using Maple 2022. 
# 
# based on original Kovacic paper description of the algorithm. 
# 
# This file contains two modules. The first is called kovacic_class 
# used to create kovacic object. The second is kovacic_tester module 
# used for unit testing the kovacic_class module 
# 
# To use this file just do 
# 
#    read "KOV.mpl" 
#    ode := diff(y(x),x$2)=2/x^2*y(x) 
#    o   := Object(kovacic_class,ode,y(x)); 
#    sol := o:-dsolve(); 
# 
# Make sure to set the currentdir() in Maple correctly in order to find 
# where you downloaded the file "KOV.mpl" to. 
#-------------------------------------------------------------------- 
 
#-------------------------------------------------------------------- 
# This class is used to solve an ode using Kovacic algorithm. 
# 
#  Please see documantion section in the paper for additional 
#  information on using this class. 
#-------------------------------------------------------------------- 
 
kovacic_class :=module() 
option object; 
 
#class to hold one entry in the gamma set for case 1 
local case_one_gamma_entry := module() 
    option object; 
    export pole_location := 0; 
    export pole_order    := 0; 
    export sqrt_r        := 0; 
    export alpha_plus    := 0; 
    export alpha_minus   := 0; 
    export b:=0; 
end module; 
 
#class to hold O_infinity information for case 1 
local case_one_O_inf := module() 
    option object; 
    export sqrt_r_inf      := 0; 
    export alpha_plus_inf  := 0; 
    export alpha_minus_inf := 0; 
    export a := 0; 
    export b := 0; 
end module; 
 
#class to hold one entry in the gamma set for case 2,3 
local case_2_and_3_gamma_entry:=module() 
   option object; 
 
   export pole_location := 0; 
   export pole_order    := 0; 
   export Ec::set       := {}; 
   export b := 0; 
end module; 
 
#-------------------------------------------------------------------- 
#  PRIVATE variables for the kovacic class 
#  only methods inside this class can access these 
#-------------------------------------------------------------------- 
 
local original_ode; 
 
#coefficients of original linear ode A y''+ B y' + C y =0 
local A,B,C; 
 
#original ode dependent and independent variables 
local y::symbol,x::symbol; 
 
local modified_ode;  #this is the z''=r*z    ode; 
local z::symbol; #dependent variable for the r_ode z''(x)=r*z(x) 
local r,s,t;     #where r=s/t; 
 
local O_inf; #order of r at infinity. degree(s)-degree(t) 
 
# poles or r. format is [ [pole,order],[pole,order],...] 
# if no poles, for example r=x, then list is empty [] 
# this means pole order zero is not in the list, since no pole. 
local poles_list::list := []; 
 
#contains all possible cases detected. Hence [1] or [1,2] or [1,2,3] 
#or remains empty [] for case 4. 
local list_of_possible_cases::list := []; 
 
local case_used_to_solve::integer:=-1; #this will have case used: 1,2 or 3 
 
#this will have n=4,6 or 12 used for case 3 only 
local n_used_for_case_3::integer:=-1; 
 
#-------------------------------------------------------------------- 
#  CONSTRUCTOR 
# 
#  the input is  the ode itself as first argument, and the 
#  dependent variable, y(x) for example, as second argument. 
#-------------------------------------------------------------------- 
export ModuleCopy::static:= proc( _self, 
            proto::kovacic_class, 
            ode::`=`,func::function(name) ,$) 
    local A,B,C; 
    local x,y,r,z::nothing,s,t; 
 
    x,y,A,B,C := parse_and_validate_ode(ode,func); 
 
    _self:-original_ode := ode; 
 
    #force r to be relatively prime ratio of 2 polynomials 
    r := normal( (2*diff(B,x)*A-2*B*diff(A,x)+B^2-4*A*C)/(4*A^2)); 
    if not type(r,'ratpoly'(anything,x)) then 
        ERROR("r= ", r ," is not polynomial or rational function in ",x); 
    fi; 
 
    s := numer(r); 
    t := denom(r); 
 
    #save all findings into the object private data 
    #so it can be used later 
    _self:-modified_ode := diff(z(x),x$2) = r*z(x); 
    _self:-O_inf := degree(t,x)-degree(s,x); 
 
    _self:-r  := r; 
    _self:-x  := x; 
    _self:-y  := y; 
    _self:-z  := z; 
    _self:-s  := s; 
    _self:-t  := t; 
    _self:-C  := C; 
    _self:-B  := B; 
    _self:-A  := A; 
 
    #determine all poles and order 
    _self:-generate_poles_and_order_list(); 
 
    #determine all possible kovacic cases 
    _self:-find_possible_cases(); 
 
    return NULL; 
end proc; 
 
#-------------------------------------------------------------------- 
# This module private function is called by constructor to parse and 
# validate the ode 
#-------------------------------------------------------------------- 
local parse_and_validate_ode:=proc(ode::`=`,func::function(name),$) 
    local x,y,A,B,C,L::list; 
    local item,dep_variables_found; 
 
    if nops(func)<>1 then 
        ERROR("dependent variable ",func," has more than one argument"); 
    fi; 
 
    y := op(0,func); 
    x := op(1,func); 
 
    #basic verification 
    if not has(ode,y) then ERROR("Supplied ",ode," has no ",y); fi; 
    if not has(ode,x) then ERROR("Supplied ", ode," has no ",x); fi; 
    if not has(ode,func) then ERROR("Supplied ",ode," has no ",func); fi; 
 
    #check it is second order 
    if PDEtools:-difforder(ode)<>2 then 
        ERROR("Only second order ode's can be used in Kocacic algorithm. "); 
    fi; 
 
    #check all dependent variables in ode which is y, match given func y(x) 
    try 
        dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); 
    catch: 
        ERROR(lastexception); 
    end try; 
 
    #o over dep_variables_found and check the 
    #independent variable is same as x i.e. ode can be y'(z)+y(z)=0 but 
    #function is y(x). 
    for item in dep_variables_found do 
        if not type(item,function) then 
            ERROR("Parse error. Expected ",func," found ",item," in", ode); 
        else 
            if op(1,item) <> x then 
                ERROR("Parse error. Expected ",func," found ",item," in",ode); 
            fi; 
            if nops(item)<>1 then 
                ERROR("Parse error. One argument allowed in ", 
                      func," found ", item," in " , ode); 
            fi; 
        fi; 
    od; 
 
    #now go over all indents in ode and check that y shows as y(x) 
    #and not as just y as the PDEtools:-Library:-GetDepVars([_self:-y],ode) 
    #code above does not detect this. i.e. it does not check  y'(x)+y=0 
    if numelems(indets(ode,identical(y))) > 0 then 
       ERROR("Parsing error, Can not have ",y," with no argument inside ",ode); 
    fi; 
 
    #check ode is linear in y(x) 
    if not has(DEtools:-odeadvisor(ode,func,['linear']),'_linear') then 
        ERROR("Only linear ode's can be used in Kocacic algorithm. "); 
    fi; 
 
    #extract coefficients of the ode 
    L:=DEtools:-convertAlg(ode,func); #this only works on linear ode's 
 
    if L[2]<>0 then ERROR("Not homogeneous ode");    fi; 
 
    #Finished parsing the ode.  A y''+B y' + C y = 0 
    C  := L[1,1]; 
    B  := L[1,2]; 
    A  := L[1,3]; 
 
    if not type(A,'ratpoly'(anything,x)) or  not type(B,'ratpoly'(anything,x)) 
        or not type(C,'ratpoly'(anything,x)) then 
        error "ode coefficients are not rational functions of ",x; 
    fi; 
 
    return x,y,A,B,C; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# main API.  Called to solve the ode using Kovacic algorithm. 
# returns the solution in the form of y(x)=.... or if no solution 
# is found returns FAIL. 
# see user guide how to use. 
#-------------------------------------------------------------------- 
export dsolve::static:=proc(_self,$) 
    local sol:=FAIL, current_case_number::posint; 
 
    if nops(_self:-list_of_possible_cases) = 0 then 
        return FAIL; 
    fi; 
 
    #keep trying all possible cases, starting from case 1 to case 3. 
    #until one case succeed or all are tried 
 
    for current_case_number in _self:-list_of_possible_cases do 
        sol := _self:-dsolve_case(current_case_number); 
        if sol<>FAIL then 
            _self:-case_used_to_solve:=current_case_number; 
            return sol; 
        fi; 
    od; 
 
    return FAIL; 
 
end proc; 
 
#-------------------------------------------------------------------- 
#  Called to solve the ode using a specific case number 
#  made public to allow user to solve using specific case 
#-------------------------------------------------------------------- 
export dsolve_case::static:=proc(_self,case_number::posint,$) 
    local sol; 
 
    if case_number>3 then ERROR("Only case number 1,2,3 are allowed");    fi; 
 
    if nops(_self:-list_of_possible_cases)=0 then 
        ERROR("No possible cases detected for this ode"); 
    fi; 
 
    if not member(case_number,_self:-list_of_possible_cases) then 
        ERROR("Case ", case_number, 
           " not one of possible cases ", _self:-list_of_possible_cases); 
    fi; 
 
    if case_number = 1 then 
       sol:= _self:-solve_case_1(); 
    elif case_number = 2 then 
       sol:= _self:-solve_case_2(); 
    else 
       sol:= _self:-solve_case_3(); 
    fi; 
 
    _self:-case_used_to_solve := case_number; 
    return sol; 
 
end proc; 
#-------------------------------------------------------------------- 
# returns back the z''(x) = r(x) z(x) ode, which is the one 
# actually solved by Kovacic algorithm 
#-------------------------------------------------------------------- 
export get_z_ode::static:=proc(_self,$) 
    local x := _self:-x; 
    local r := _self:-r; 
    local z := _self:-z; 
 
    if _self:-s = 0 then 
         return diff(z(x),x$2) = 0; 
    else 
         return diff(z(x),x$2) = numer(r)%/denom(r) * z(x); 
    fi; 
end proc; 
 
#-------------------------------------------------------------------- 
# returns back the r term in the z''(x) = r(x) z(x) ode 
#-------------------------------------------------------------------- 
export get_r::static:=proc(_self,$) 
    if _self:-s = 0 then 
        return 0; 
    else 
        return numer(_self:-r)%/denom(_self:-r); 
    fi; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns s, where r = s/t and z'' = r z(t) 
#-------------------------------------------------------------------- 
export get_s::static:=proc(_self,$) 
    return _self:-s; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns t, where r = s/t and z'' = r z(t) 
#-------------------------------------------------------------------- 
export get_t::static:=proc(_self,$) 
    return _self:-t; 
end proc; 
 
#-------------------------------------------------------------------- 
# Returns back the original ode ( A y''+ B y' + C y = 0 ) 
#-------------------------------------------------------------------- 
export get_y_ode::static:=proc(_self,$) 
    return _self:-original_ode; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns list of poles of r, where  z''(x) = r z(x) 
#  The list has format [ [pole location,order], ...] 
#-------------------------------------------------------------------- 
export get_poles::static:=proc(_self,$) 
    return _self:-poles_list; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns O_infinity of r, where  z''(x) = r z(x) 
#-------------------------------------------------------------------- 
export get_order_at_infinity::static:=proc(_self,$) 
    return _self:-O_inf; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns list of possible kovacic cases possible. 
#-------------------------------------------------------------------- 
export get_possible_cases::static:=proc(_self,$)::list; 
    return _self:-list_of_possible_cases; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns actual case number used when solving ode. 
#  can be 1,2 or 3 
#  if no cases applicable, then -1 is returned 
#-------------------------------------------------------------------- 
export get_case_used::static:=proc(_self,$)::integer; 
    return _self:-case_used_to_solve; 
end proc; 
 
#-------------------------------------------------------------------- 
#  returns n used for case 3. Either 4,6, or 12. 
#  if not case 3 used, then -1 is returned. 
#-------------------------------------------------------------------- 
export get_n_case_3::static:=proc(_self,$)::integer; 
    return _self:-n_used_for_case_3; 
end proc; 
 
 
#-------------------------------------------------------------------- 
# All functions below are private 
#-------------------------------------------------------------------- 
 
#-------------------------------------------------------------------- 
#This proc find all possible kovacic cases. These can be 1,2 or 3. 
#if none of these found, then empty list is returned, which is 
#case 4 in the paper. For example, if case 1 is only possible, 
#then [1] is returned. if case 1 and 2 are possible, then [1,2] 
#is returned, if all three cases possible, then [1,2,3] is returned. 
#If no cases possible then [] returned. 
#-------------------------------------------------------------------- 
local find_possible_cases::static:=proc(_self,$) 
    local L::list := []; 
    local poles_order := convert(_self:-poles_list[..,2],set); 
    local T::set; 
 
    #check for case 1 
    T := select(Z-> Z>2,poles_order); 
    if  nops( select(Z->type(Z,odd),T) )=0  and 
        (_self:-O_inf<0 and type(_self:-O_inf,even)) or 
            _self:-O_inf=0 or _self:-O_inf>1 then 
        L:= [1]; 
    fi; 
 
    if nops(poles_order)>0 then #must have at least one pole for 2,3 cases 
 
        if nops(poles_order)<>0 then #can not have pole order 0, i.e. no poles 
 
            T:=select(Z-> Z>2,poles_order); 
 
            # r must have at least one pole that is either odd order greater 
            #than 2 or else has order 2 
 
            if  nops( select(Z->type(Z,odd),T) )<>0 
                     or member(2,poles_order) then 
                L:= [ op(L),2 ]; 
            fi; 
 
            #check for case 3 
 
            if not member(0,poles_order) 
                   and nops(select(Z-> Z>2,poles_order))=0 
                   and _self:-O_inf>=2 then 
 
               L:= [ op(L),3 ]; 
 
            fi; 
        fi; 
    fi; 
 
    _self:-list_of_possible_cases := L; 
 
    return NULL; 
end proc; 
 
#-------------------------------------------------------------------- 
#called to find all poles of r and the order 
#of each pole. Uses Maple's sqrfree. 
#-------------------------------------------------------------------- 
local generate_poles_and_order_list::static:=proc(_self,$) 
    local L::list := []; 
    local sol::list; 
    local current_sol; 
    local poles_list::list; 
    local current_pole; 
    local r := _self:-r, x := _self:-x; 
 
    poles_list := sqrfree(denom(r),x); 
    poles_list := poles_list[2,..]; #we do not need the overall factor 
 
    if nops(poles_list) = 0 then 
        _self:-poles_list:=[]; 
    else 
        for current_pole in poles_list do 
 
            sol := solve(current_pole[1]=0,[x]); 
            sol := ListTools:-Flatten(sol); 
 
            for current_sol in sol do 
                L := [op(L), [rhs(current_sol) ,current_pole[2] ] ]; 
            od; 
 
        od; 
 
        _self:-poles_list := L; 
    fi; 
 
    return NULL; 
 
end proc: 
#-------------------------------------------------------------------- 
# 
#       C A S E     O N E      I M P L E M E N T A T I O N 
# 
#  returns ode solution using case 1, or FAIL is no solution exist 
#-------------------------------------------------------------------- 
local solve_case_1::static:=proc(_self,$) 
    local O_infinity_set::kovacic_class:-case_one_O_inf; 
    local gamma_set::set(kovacic_class:-case_one_gamma_entry):={}; 
 
    gamma_set, O_infinity_set := _self:-case_1_step_1(); 
    return _self:-case_1_step_2(gamma_set, O_infinity_set); 
end proc; 
 
#-------------------------------------------------------------------- 
# called from _self:-solve_case_1() 
#-------------------------------------------------------------------- 
local case_1_step_1::static:=proc(_self,$):: 
          set(kovacic_class:-case_one_gamma_entry), 
          kovacic_class:-case_one_O_inf; 
 
    local current_pole; 
    local e::kovacic_class:-case_one_gamma_entry; 
    local o::kovacic_class:-case_one_O_inf; 
    local b,a,v,i::integer; 
    local b_coeff_in_r,b_coeff_in_r_inf_square; 
    local x := _self:-x, r := _self:-r; 
    local N::integer; 
    local laurent_c; 
    local current_term; 
    local b_from_r,b_from_laurent_series; 
 
    #this contains all information generated for each pole of r 
    #this is what is called the set GAMMA in the paper and 
    #in the diagram of algorithm above. 
    local gamma_set::set(kovacic_class:-case_one_gamma_entry) := {}; 
 
    if nops(_self:-poles_list) = 0 then 
        gamma_set := {}; 
    else 
        for current_pole in _self:-poles_list do 
            e := Object(kovacic_class:-case_one_gamma_entry); 
 
            e:-pole_location  := current_pole[1]; 
            e:-pole_order     := current_pole[2]; 
 
            if e:-pole_order =1 then 
 
               e:-sqrt_r      := 0; 
               e:-alpha_plus  := 1; 
               e:-alpha_minus := 1; 
               gamma_set      := { op(gamma_set), e }; 
 
            elif e:-pole_order = 2 then 
 
                e:-sqrt_r := 0; 
                e:-b      := _self:-b_partial_fraction(r,x,e:-pole_location,2); 
                e:-alpha_plus  := 1/2+1/2*sqrt(1+4*e:-b); 
                e:-alpha_minus := 1/2-1/2*sqrt(1+4*e:-b);; 
                gamma_set := { op(gamma_set), e }; 
 
            else 
 
                v         := e:-pole_order/2; 
                e:-sqrt_r := 0; 
 
                for N from 2 to v  do 
                    laurent_c := _self:-laurent_coeff(sqrt(r),x, 
                                                      e:-pole_location,v,N); 
                    current_term := laurent_c/(x-e:-pole_location)^N; 
                    e:-sqrt_r    := e:-sqrt_r + current_term; 
                    if N = v then 
                        a := laurent_c; 
                    fi; 
                od; 
 
                b_from_r :=_self:-b_partial_fraction(r,x,e:-pole_location,v+1); 
                b_from_laurent_series :=_self:-laurent_coeff(sqrt(r),x, 
                                                    e:-pole_location,v,v+1); 
 
                e:-b := b_from_r - b_from_laurent_series; 
 
                e:-alpha_plus  := 1/2*((e:-b)/a  + v); 
                e:-alpha_minus := 1/2*(-(e:-b)/a + v); 
                gamma_set      := { op(gamma_set), e }; 
 
            fi; 
        od; 
    fi; 
 
    o := Object(case_one_O_inf); 
 
    if _self:-O_inf > 2 then 
 
        o:-sqrt_r_inf      := 0; 
        o:-alpha_plus_inf  := 0; 
        o:-alpha_minus_inf := 1; 
 
    elif _self:-O_inf=2 then 
 
        o:-sqrt_r_inf :=0; 
        o:-b          := lcoeff(_self:-s) / lcoeff(_self:-t); 
        b             := radsimp((1+4*o:-b)^(1/2)); 
        o:-alpha_plus_inf  := 1/2 + 1/2*b; 
        o:-alpha_minus_inf := 1/2 - 1/2*b; 
 
    else #order at infinity -2*v<= 0 which must be even 
 
        v     := (-_self:-O_inf) / 2; 
        o:-sqrt_r_inf :=0; 
 
        for i from 0 to v do 
 
            laurent_c     := _self:-laurent_coeff(sqrt(r),x,infinity,v,i); 
            o:-sqrt_r_inf := o:-sqrt_r_inf + laurent_c*x^i; 
            if i = v then 
               o:-a := laurent_c; 
            fi; 
 
        od; 
 
        b_coeff_in_r_inf_square := _self:-laurent_coeff( 
                                            o:-sqrt_r_inf^2,x,0,v,v-1); 
 
        b_coeff_in_r       := _self:-get_coefficient_of_r(v-1); 
        o:-b               := b_coeff_in_r - b_coeff_in_r_inf_square; 
        o:-alpha_plus_inf  := 1/2*(  (o:-b)/(o:-a) - v); 
        o:-alpha_minus_inf := 1/2*( -(o:-b)/(o:-a) - v); 
    fi; 
 
    return gamma_set,o; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# Finds b coefficient in r for the case one only when v<=0 
# using long division 
#-------------------------------------------------------------------- 
local get_coefficient_of_r::static:=proc(_self,the_degree,$) 
    local r := _self:-r; 
    local x := _self:-x; 
    local c,R,t; 
 
    try 
        c := coeff(r,x,the_degree); 
    catch: 
        R := rem(numer(r),denom(r),x); 
        t := denom(r); 
        c := lcoeff(R)/lcoeff(t); 
    end try; 
 
    return c; 
end proc; 
 
#-------------------------------------------------------------------- 
# called from _self:-solve_case_1() 
# This determines the set of d non-negative integers and 
# corresponding w for each d. 
#-------------------------------------------------------------------- 
local case_1_step_2::static:=proc(_self, 
           gamma_set::set(kovacic_class:-case_one_gamma_entry), 
           O_infinity::kovacic_class:-case_one_O_inf,$) 
 
    local d,N,K,number_of_poles::integer,current_alpha_infinity; 
    local item::list; 
    local the_sign::integer,sign_list::list; 
    local w; 
    local r_solution,y_solution; 
    local x := _self:-x; 
 
    #this will contains all data found for any nonnegative d. Each 
    #entry will be a list of this form 
    #       [d, w] 
 
    local good_d_found := Array(1..0); 
    local B::Matrix; #this will contain good_d_found but as matrix 
 
    number_of_poles := nops(gamma_set); 
    sign_list := combinat:-permute([1$number_of_poles,-1$number_of_poles], 
                                   number_of_poles); 
 
    for K,current_alpha_infinity in 
                [O_infinity:-alpha_plus_inf,O_infinity:-alpha_minus_inf] do 
 
        for item in sign_list do 
 
            d := 0; 
 
            if number_of_poles>0 then 
 
                for N,the_sign in item do 
                    if the_sign = -1 then 
                        d := d-gamma_set[N]:-alpha_minus; 
                    else 
                        d := d-gamma_set[N]:-alpha_plus; 
                    fi; 
                od; 
 
            fi; 
 
            d := simplify(d + current_alpha_infinity); 
 
            if type(d,'integer') and d >= 0 then 
                w := 0; 
 
                for N,the_sign in item do 
 
                    if number_of_poles>0 then 
 
                        if the_sign = -1 then 
                            w := w +(-1)*gamma_set[N]:-sqrt_r + 
                                    (gamma_set[N]:-alpha_minus)/ 
                                          (x - gamma_set[N]:-pole_location); 
                        else 
                            w := w + gamma_set[N]:-sqrt_r  + 
                                    (  gamma_set[N]:-alpha_plus)/ 
                                         (x  - gamma_set[N]:-pole_location); 
                        fi; 
 
                    fi; 
 
                od; 
 
                #this to get the sign in according to paper. First + then - 
                if K = 1 then 
                   w := w + O_infinity:-sqrt_r_inf; 
                else 
                   w := w - O_infinity:-sqrt_r_inf; 
                fi; 
 
                good_d_found ,= [ d, w]; 
            fi; 
        od; 
    od; 
 
    if numelems(good_d_found) = 0 then return FAIL; fi; 
 
    #now convert the array to Matrix, and sort on d, so we 
    #start will the smallest 
    #degree d, which is the first column, as that will be most efficient. 
    #convert to set first, to remove any possible duplicat entries 
    #then convert to Matrix 
 
    convert(good_d_found,set); 
    B := convert(convert(%,list),Matrix); 
    B := B[sort(B[.., 1], 'output'= 'permutation')]; 
 
    for N from 1 to LinearAlgebra:-RowDimension(B) do 
 
        r_solution := _self:-case_1_step_3(B[N,1], B[N,2]); #(d,w) 
 
        if r_solution <> FAIL then 
            y_solution := _self:-build_y_solution_from_r_solution(r_solution); 
            return y_solution; 
        fi; 
 
    od; 
 
    return FAIL; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# called from _self:-case_1_step_2() 
#-------------------------------------------------------------------- 
local case_1_step_3::static:=proc(_self,d::nonnegative,w,$) 
    local x := _self:-x; 
    local r := _self:-r; 
    local p; 
    local i::integer; 
    local a::nothing; 
    local eq; 
    local coeff_sol; 
    local tmp,W; 
    local final_result; 
 
    p := x^d; 
 
    for i from d-1 by -1 to 0 do 
        p := p + a[i] * x^i; 
    od; 
 
    #using original kovacis method. Not Smith. 
    eq := simplify( diff(p,x$2)+2*w*diff(p,x)+(diff(w,x)+w^2-r)*p) = 0; 
 
    if d = 0 then 
        if not evalb(eq) then return FAIL; fi; 
    else 
        #solve for coefficients 
        try 
           coeff_sol := timelimit(30,solve( 
                                   identity(eq,x), [seq(a[i],i=0..d-1)])); 
        catch: 
            return FAIL; 
        end try; 
 
        if nops(coeff_sol) = 0 then return FAIL; fi; 
 
        tmp := map(evalb,coeff_sol[1]); 
        if has(tmp,true) then  return FAIL; fi; 
 
        p := eval(p,coeff_sol[1]); #to force a[i] solutions to update 
    fi; 
 
    W := diff(p,x)/p + w; 
    W := radsimp(W); 
 
    tmp := diff(W,x)+W^2; 
 
    if evalb( tmp = r) or is(tmp = r) or simplify(tmp-r)=0 then #can be used 
        try 
            tmp := int(w, x); 
        catch: 
            return FAIL; 
        end try; 
 
        if has(tmp,int) then return FAIL; fi; 
 
        final_result := simplify(p*exp(tmp)); 
        if has(final_result,signum) or has(final_result,csgn) then 
            final_result := p*exp(tmp); 
        fi; 
 
        return final_result; 
    else 
        return FAIL; 
    fi; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# 
#  C A S E    T W O   I M P L E M E N T A T I O N 
# 
#  returns ode solution using case 2, or FAIL is no solution exist 
#-------------------------------------------------------------------- 
local solve_case_2::static:=proc(_self,$) 
 
    local E_inf::set; 
    local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry):={}; 
 
    gamma_set, E_inf := _self:-case_2_step_1(); 
    return _self:-case_2_step_2(gamma_set, E_inf ); 
 
end proc; 
 
#-------------------------------------------------------------------- 
# 
#-------------------------------------------------------------------- 
local case_2_step_1::static:=proc(_self,$):: 
        set(kovacic_class:-case_2_and_3_gamma_entry),set; 
 
    local current_pole; 
    local e::kovacic_class:-case_2_and_3_gamma_entry; 
    local E_inf::set; 
    local b; 
    local x := _self:-x; 
    local r := _self:-r; 
 
    #this contains all information generated for each pole of r 
    #this is what is called the set GAMMA in the paper and in the diagram of 
    #algorithm above. 
    local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry) := {}; 
 
    for current_pole in _self:-poles_list do 
 
        e := Object(kovacic_class:-case_2_and_3_gamma_entry); 
        e:-pole_location  := current_pole[1]; 
        e:-pole_order     := current_pole[2]; 
 
        if e:-pole_order = 1 then 
 
           e:-Ec     := {4}; 
           gamma_set := { op(gamma_set), e }; 
 
        elif e:-pole_order = 2 then 
 
            e:-b      := _self:-b_partial_fraction(r,x,e:-pole_location,2); 
            e:-Ec     := {2,2+2*sqrt(1+4* e:-b),2-2*sqrt(1+4* e:-b)}; 
            e:-Ec     := select(z->type(z,integer),e:-Ec); 
            gamma_set := { op(gamma_set), e }; 
 
        else 
 
           e:-Ec     := {e:-pole_order}; 
           gamma_set := { op(gamma_set), e }; 
 
        fi; 
    od; 
 
    if _self:-O_inf>2 then 
 
        E_inf := {0,2,4}; 
 
    elif _self:-O_inf=2 then 
 
        b     := lcoeff(_self:-s) / lcoeff(_self:-t); 
        E_inf := {2,2+2*sqrt(1+4*b),2-2*sqrt(1+4*b)}; 
        E_inf := select(z->type(z,integer),E_inf); 
 
    else #order at infinity v< 2 
 
        E_inf := {_self:-O_inf}; 
 
    fi; 
 
    return gamma_set,E_inf; 
end proc; 
 
#-------------------------------------------------------------------- 
# called from _self:-solve_case_2() 
# This determines the set of d non-negative integers and 
# corresponding w for each d. 
#-------------------------------------------------------------------- 
local case_2_step_2::static:=proc(_self, 
           gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry), 
           E_inf::set, 
           $) 
 
    local L::list := []; 
    local item,current_E_inf; 
    local d,N,ee,theta; 
    local x := _self:-x; 
    local r_solution,y_solution; 
 
    #this will contains all data found for any nonnegative d. Each 
    #entry will be a list of this form 
    #       [d, theta] 
    #so if we obtain say 3 values of d that are nonnegative, 
    #there will be 3 such lists in this array 
 
    local good_d_found := Array(1..0); 
    local B::Matrix; #this will contain good_d_found as matrix 
 
    for item in gamma_set do 
        L := [ op(L), convert(item:-Ec,list) ]; 
    od; 
 
    #now find all possible tuples 
    if nops(L)>1 then 
       L := kovacic_class:-cartProdSeq(op(L)); 
    else 
       L := L[1]; 
    fi; 
 
    #DEBUG(); 
    for current_E_inf in E_inf do 
        for item in L do 
 
            d := 1/2*( current_E_inf - add(item)); 
 
            if type(d,'integer') and d >= 0 then 
                theta :=0; 
 
                for N,ee in item do 
                    theta := theta + ee/(x-gamma_set[N]:-pole_location); 
                od; 
 
                theta         := 1/2*theta; 
                good_d_found ,= [ d, theta]; 
            fi; 
 
        od; 
    od; 
 
    if numelems(good_d_found) = 0 then return FAIL; fi; 
 
    #now convert the array to Matrix, and sort on d, so 
    #we start will the smallest degree d, which is the first column, as 
    #that will be most efficient. 
    #convert to set first, to remove any possible duplicat entries 
    #then convert to Matrix 
 
    convert(good_d_found,set); 
    B := convert(convert(%,list),Matrix); 
    B := B[sort(B[.., 1], 'output'= 'permutation')]; 
 
    for N from 1 to LinearAlgebra:-RowDimension(B) do 
 
        r_solution := _self:-case_2_step_3(B[N,1], B[N,2]);  #(d,theta) 
 
        if r_solution <> FAIL then 
            y_solution := _self:-build_y_solution_from_r_solution(r_solution); 
            return y_solution; 
        fi; 
 
    od; 
 
    return FAIL; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# called from _self:-case_2_step_2() 
#-------------------------------------------------------------------- 
local case_2_step_3::static:=proc(_self,d::nonnegative,theta,$) 
 
    local p, i; 
    local a::nothing; 
    local tmp; 
    local coeff_sol := []; 
    local eq; 
    local phi; 
    local sol_w; 
    local current_w; 
    local r:=_self:-r; 
    local w::nothing; 
    local x := _self:-x; 
    local sol; 
 
    p := x^d; 
 
    for i from d-1 by -1 to 0 do 
        p := p + a[i] * x^i; 
    od; 
 
    eq:= simplify(diff(p,x$3) + 3*theta*diff(p,x$2) + 
            (3*theta^2+3*diff(theta,x) -4*r) * diff(p,x) 
              + ( diff(theta,x$2)+3*theta*diff(theta,x)+theta^3 - 
                 4*r*theta - 2*diff(r,x))*p) = 0; 
 
    if d=0 then 
        if not evalb(eq) then return FAIL; fi; 
    else 
        #solve for polynomial coefficients 
        try 
            coeff_sol:= timelimit(30,solve( 
                                   identity(eq,x), [seq(a[i],i=0..d-1)])); 
        catch: 
            return FAIL; 
        end try; 
 
        if nops(coeff_sol) = 0 then return FAIL; fi; 
 
        tmp := map(evalb,coeff_sol[1]); 
        if has(tmp,true) then return FAIL; fi; 
 
        p := eval(p,coeff_sol[1]); #to force a[i] solutions to update 
    fi; 
 
    phi := theta + diff(p,x)/p; 
    eq  := w^2 - phi*w +  simplify(1/2*diff(phi,x)+1/2*phi^2-r) = 0; 
 
    try 
        sol_w := timelimit(30,solve(eq,[w])); #changed to [] 
        sol_w := ListTools:-Flatten(sol_w); 
    catch: 
        return FAIL; 
    end try; 
 
    if nops(sol_w) = 0 then return FAIL; fi; 
 
    for current_w in sol_w do 
 
        current_w := radsimp(rhs(current_w)); 
 
        #verify w before using it. Added 1/12/2022 4 PM 
        tmp := diff(current_w,x)+current_w^2; 
        if evalb( tmp= r) or is(tmp = r) or simplify(tmp-r)=0 then 
            try 
               sol := timelimit(30,int(current_w, x)); 
            catch: 
               return FAIL; 
            end try; 
 
            if has(sol,int) then return FAIL; fi; 
 
            return simplify(exp(sol)); 
        fi; 
 
    od; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# 
#  C A S E   T H R E E    I M P L E M E N T A T I O N 
# 
#  returns ode solution using case 3, or FAIL is no solution exist 
#-------------------------------------------------------------------- 
 
local solve_case_3::static:=proc(_self,$) 
 
    local E_inf::set; 
    local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry):={}; 
    local sol; 
 
    #these are possible degress of w to try until one works or none works 
    #local w_degree::list := [4,6,12]; 
    local w_degree::list := [4,6,12]; 
    local n::posint; 
 
    for n in w_degree do 
        gamma_set, E_inf := _self:-case_3_step_1(n); 
        sol              := _self:-case_3_step_2(gamma_set, E_inf, n ); 
 
        if sol<>FAIL then return sol; fi; 
    od; 
 
    return FAIL; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# First step in case 3 
#-------------------------------------------------------------------- 
local case_3_step_1::static:=proc(_self, 
           n::posint, 
           $)::set(kovacic_class:-case_2_and_3_gamma_entry),set; 
 
 
    local current_pole; 
    local e::kovacic_class:-case_2_and_3_gamma_entry; 
    local E_inf::set; 
    local b,k; 
    local x := _self:-x; 
    local r := _self:-r; 
 
    #this contains all information generated for each pole of r 
    #this is what is called the set GAMMA in the paper and in the diagram of 
    #algorithm above. 
    local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry):={}; 
 
    for current_pole in _self:-poles_list do 
 
        e := Object(kovacic_class:-case_2_and_3_gamma_entry); 
 
        e:-pole_location := current_pole[1]; 
        e:-pole_order    := current_pole[2]; 
 
        if e:-pole_order = 1 then 
 
           e:-Ec     := {12}; 
           gamma_set := { op(gamma_set), e }; 
 
        elif e:-pole_order = 2 then 
 
            e:-b      := _self:-b_partial_fraction(r,x,e:-pole_location,2); 
            e:-Ec     := {seq( 6+12*k/n*sqrt(1+4*(e:-b)),k=-n/2..n/2,1)}; 
            e:-Ec     := select(z->type(z,integer),e:-Ec); 
            gamma_set := { op(gamma_set), e }; 
 
        else 
           ERROR("Internal error. Case 3 can only have poles of order 1 or 2"); 
        fi; 
    od; 
 
    #same formula, but different b 
    b     := lcoeff(_self:-s) / lcoeff(_self:-t); 
    E_inf := {seq( 6+12*k/n*sqrt(1+4*b),k=-n/2..n/2,1)}; 
    E_inf := select(z->type(z,integer),E_inf); 
 
    return gamma_set,E_inf; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# Second step in case 3 
# 
# called from _self:-solve_case_3() 
# This determines the set of d non-negative integers and 
# corresponding w for each d. 
#-------------------------------------------------------------------- 
local case_3_step_2::static:=proc( 
        _self, 
        gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry), 
        E_inf::set, 
        n::posint,$) 
 
    local L::list := []; 
    local item,current_E_inf; 
    local d,ee,theta,N; 
    local x := _self:-x; 
    local r_solution,y_solution; 
    local S; 
    local current_iteration::integer; 
    local tmp; 
 
    #this will contains all data found for any nonnegative d. Each 
    #entry will be a list of this form 
    #       [d, theta, S] 
    #so if we obtain say 3 values of d that are nonnegative 
    local good_d_found := Array(1..0); 
    local B::Matrix; #this will contain good_d_found as matrix 
 
    #DEBUG(); 
    for item in gamma_set do 
        L := [ op(L), convert(item:-Ec,list) ]; 
    od; 
 
    #now find all possible tuples 
    if nops(L)>1 then 
       L := kovacic_class:-cartProdSeq(op(L)); 
    else 
       L := L[1]; 
    fi; 
 
    current_iteration:=0; 
    for current_E_inf in E_inf do 
        for item in L do 
            current_iteration := current_iteration+1; 
 
            d := n/12*( current_E_inf - add(item)); 
 
            if type(d,'integer') and d >= 0 then 
                theta :=0; 
                for N,ee in item do 
                    theta := theta + ee/(x-gamma_set[N]:-pole_location); 
                od; 
                theta := n/12*theta; 
                theta := simplify(theta); 
                S := mul( (x-gamma_set[N]:-pole_location), N=1..nops(item)); 
                tmp := simplify(S) assuming real; 
                if not has(tmp,csgn) and not has(tmp,signum) then 
                   S:= tmp; 
                fi; 
 
                good_d_found ,= [ d, theta, S]; 
            fi; 
        od; 
    od; 
 
    #now convert the array to Matrix, and sort on d, so we start 
    # with the smallest 
    #degree d, which is the first column, as that will be most efficient. 
 
    if numelems(good_d_found) = 0 then return FAIL; fi; 
 
    #convert to set first, to remove any possible duplicat entries 
    #then convert to Matrix 
    convert(good_d_found,set); 
    B := convert(convert(%,list),Matrix); 
    B := B[sort(B[.., 1], 'output'= 'permutation')]; 
 
    for N from 1 to LinearAlgebra:-RowDimension(B) do 
        #(d,theta,S,n) 
        r_solution := _self:-case_3_step_3(B[N,1], B[N,2], B[N,3],n); 
 
        if r_solution <> FAIL then 
            _self:-n_used_for_case_3 := n; 
            y_solution := _self:-build_y_solution_from_r_solution(r_solution); 
            return y_solution; 
        fi; 
    od; 
 
    return FAIL; 
end proc; 
 
#-------------------------------------------------------------------- 
# Third and final step in case 3. 
# called from _self:-case_3_step_2(). Returns solution or FAIL is no 
# solution found. 
#-------------------------------------------------------------------- 
local case_3_step_3::static:=proc( 
                 _self, 
                 d::nonnegative, 
                 theta, 
                 S, 
                 n::posint, 
                 $) 
 
    local p, i,P_minus_1; 
    local a::nothing; 
    local final_result,tmp; 
    local coeff_sol := []; 
    local sol_w; 
    local current_w; 
    local r := _self:-r; 
    local x := _self:-x; 
    local omega::symbol; 
    local P := Array(-1..n); 
    local result_of_simplify; 
    local result_of_is_check; 
    local omega_equation; 
 
    #this makes p(x). For example if d=3, then the result will be 
    #         p(x) = x^3 + a(2) x^2 + a(1) x  + a(0) 
    #where the number of unknowns to determine is always the same 
    #as the degree, in this case a(0),a(1),a(2) 
 
    p := x^d; #this will be 1 if degree is zero 
 
    for i from d-1 by -1 to 0 do 
        p := p + a[i] * x^i; 
    od; 
 
    #build the P_n polynomials based on p(x) above 
    P[n] := -p; 
 
    for i from n by -1 to 0 do 
 
        if i=n then 
            P[i-1] := -S * diff(P[i], x) - S*theta*P[i]; 
            P[i-1] := simplify( P[i-1]); 
        else 
            P[i-1] := -S * diff(P[i], x) + 
                  ( (n-i)*diff(S,x) - S*theta)*P[i] - (n-i)*(i+1)*S^2*r*P[i+1]; 
            P[i-1] := simplify( P[i-1]); 
        fi; 
 
    od; 
 
    #solve P[-1] for p(x) if needed. 
    P_minus_1 := expand(numer(radsimp(P[-1]))); 
 
    if P_minus_1 <>0 and evalb(p<>1) then 
        if not hastype(P_minus_1,'indexed')  then #it must be indexed now 
           ERROR("Internal error. Please report. This should not happen"); 
        fi; 
 
        try 
            coeff_sol := timelimit(30,solve(identity(P_minus_1,x), [seq(a[i],i=0..d-1)])); 
        catch: 
            return FAIL; 
        end try; 
 
        if nops(coeff_sol)=0 then #unable to solve 
           return FAIL; 
        fi; 
 
        map(evalb,coeff_sol[1]); #check all solved for 
        if has(%,true) then return FAIL; fi; 
    fi; 
 
    #build the equation for omega 
    omega_equation := 0; 
 
    for i from 0 to n do 
        omega_equation := omega_equation + S^i*P[i]/(n-i)! * omega^i ; 
    od; 
 
    omega_equation := simplify(omega_equation); 
 
    if nops(coeff_sol)<>0 then 
        #to force a[i] solutions to update to solved coefficients 
        omega_equation := eval(omega_equation,coeff_sol[1]); 
    fi; 
 
    try 
        sol_w := timelimit(30, solve(omega_equation=0, [omega])); 
        sol_w := ListTools:-Flatten(sol_w); 
    catch: 
        return FAIL; 
    end try; 
 
    if nops(sol_w) = 0 then return FAIL; fi; 
 
    #go over each w solution and use one that works. verify before using 
 
    for current_w in sol_w do 
 
        current_w := rhs(current_w); 
 
        if not has(current_w, RootOf) then 
            try 
                current_w := timelimit(30,simplify(current_w)); 
            catch: 
                NULL; 
            end try; 
 
            if is_w_verified(current_w,x,r) then 
                final_result := _self:-simplify_final_result(current_w,x); 
                if final_result<>FAIL then 
                   return final_result; 
                fi; 
            fi; 
        else #no Rootof, try to resolve 
            try 
                current_w := timelimit(30,[allvalues(current_w)]); 
            catch: 
                return FAIL; 
            end try; 
 
            current_w := current_w[1]; #just use any root. Pick first 
 
            if not has(current_w, RootOf) then 
                try 
                    current_w := timelimit(30,simplify(current_w)); 
                catch: 
                    NULL; 
                end try; 
 
                if is_w_verified(current_w,x,r) then 
                    final_result := _self:-simplify_final_result(current_w,x); 
                    if final_result<>FAIL then 
                       return final_result; 
                    fi; 
 
                fi; 
            fi; 
        fi; 
    od; 
 
    return FAIL; 
 
end proc; 
 
 
#-------------------------------------------------------------------- 
# Called from case 3, step 3 to verify w 
#-------------------------------------------------------------------- 
local is_w_verified:=proc(current_w,x,r)::truefalse; 
    local tmp; 
    local result_of_is_check::truefalse; 
    local result_of_simplify; 
 
    tmp := diff(current_w,x)+current_w^2; 
 
    if evalb( tmp=r) then return true; fi; 
 
    try 
        result_of_simplify := timelimit(30,simplify(tmp-r)); 
        if evalb(result_of_simplify=0) then return true; fi; 
    catch: 
        NULL; 
    end try; 
 
    try 
        result_of_is_check := timelimit(30,is(tmp = r)); 
        if result_of_is_check then return true; fi; 
    catch: 
        NULL; 
    end try; 
 
 
    return false; 
end proc; 
 
#-------------------------------------------------------------------- 
# Called from case 3, step 3 to simplify final result 
#-------------------------------------------------------------------- 
 
local simplify_final_result::static:=proc(_self,omega,x,$) 
    local final_result; 
    local integral_result; 
 
    try 
        integral_result := timelimit(30,int(omega, x)); 
    catch: 
        return FAIL; 
    end try; 
 
    if has(integral_result,int) or has(integral_result,Int) then 
        return exp(Int(omega, x)); 
    fi; 
 
    try 
        final_result := timelimit(30,simplify(exp(integral_result))) assuming real; 
        if has(final_result,signum) or has(final_result,csgn) 
            or has(final_result,abs) then 
            final_result := exp(integral_result); 
        fi; 
    catch: 
        final_result:= exp(integral_result); 
    end try; 
 
    return final_result; 
end proc; 
 
 
#-------------------------------------------------------------------- 
# External proc helper 
# provided thanks to Joseph Riel 
#-------------------------------------------------------------------- 
local cartProdSeq:= proc(L::seq(list)) 
    local Seq::nothing,i::nothing,j; 
    option `Copyright (C) 2007, Joseph Riel. All rights reserved.`; 
        eval([subs(Seq= seq, foldl(Seq, [cat(i, 1..nargs)], 
                                  seq(cat(i,j)= L[j], j= nargs..1, -1)))]) 
end proc: 
 
#-------------------------------------------------------------------- 
# Checks for special math function. 
# This function was provided thanks to Carl Love. 
# modified to allow some special functions. 
#-------------------------------------------------------------------- 
local has_special_math_function:= subs( 
    _F= {(op@FunctionAdvisor)~(FunctionAdvisor("class_members", "quiet"), "quiet")[]} 
        minus (   {FunctionAdvisor("elementary", "quiet")[]} 
                  union {erf,erfc, erfi, Im,Re,signum,max,argument} ), 
    proc(e::algebraic, x::{name, set(name)}:= {}) 
        hastype(e, And(specfunc(_F), dependent(x))) 
    end proc 
): 
 
#-------------------------------------------------------------------- 
# called if step 3 is successful. Build y solution from r solution 
#-------------------------------------------------------------------- 
local build_y_solution_from_r_solution::static:=proc(_self,r_solution,$) 
 
    local y1,y2; 
    local int_B_over_A; 
    local A :=_self:-A, B:=_self:-B; 
    local x :=_self:-x, y:= _self:-y; 
    local tmp; 
 
    if B = 0 then 
       y1 := r_solution; 
       int_B_over_A := 0; 
    else 
        try 
            int_B_over_A := timelimit(60,int(-B/A,x)); 
            if has(int_B_over_A,int) or has(int_B_over_A,RootOf) then 
                int_B_over_A := Int(-B/A, x); 
                y1 :=r_solution*exp(1/2*int_B_over_A); 
            else 
                y1 := simplify(r_solution*exp(1/2*int_B_over_A)); 
                if has(y1,signum) or has(y1,csgn) then 
                    y1 := r_solution*exp(1/2*int_B_over_A); 
                fi; 
            fi; 
        catch: 
            int_B_over_A := Int(-B/A,x); 
            y1 := r_solution*exp(1/2*int_B_over_A); 
        end try; 
    fi; 
 
    try 
        y2 := timelimit(30, int( simplify(exp(int_B_over_A))/y1^2,x)); 
 
        if kovacic_class:-has_special_math_function(y2,x) then 
            y2 := y1 * Int( exp(int_B_over_A)/y1^2,x) 
        else 
            y2 := y1 * y2; 
        fi; 
        tmp := simplify(y2); 
        if has(tmp,signum) or has(tmp,csgn) then 
            y2 := y1 * y2; 
        else 
            y2 := tmp; 
        fi; 
    catch: 
        y2 := y1 * Int( exp(int_B_over_A)/y1^2,x) 
    end try; 
 
    return y(x) = _C1*y1 + _C2*y2; 
 
end proc; 
 
#-------------------------------------------------------------------- 
# Helper private function called to obtain b from the partial fraction 
# decomposition of r needed by the differenent case implementations 
#-------------------------------------------------------------------- 
local b_partial_fraction::static:=proc(_self, 
        r, 
        x::symbol, 
        pole_location, 
        the_power::posint,$) 
 
    local r_partial_fraction,T::nothing; 
 
    r_partial_fraction := allvalues(convert(r,'fullparfrac',x)); 
 
    #adding dummy T so that select below always works 
    r_partial_fraction := T+r_partial_fraction; 
    select(z->hastype(z, 
            anything/(anything*identical(x - pole_location)^the_power)), 
            r_partial_fraction); 
 
    return coeff(%,1/(x-pole_location)^the_power); 
 
end proc; 
 
#-------------------------------------------------------------------- 
# Helper private function called to return Laurent series coefficient 
# of the 1/(x-c)^n  term. where n=1,2,3,.... for the function f(x) 
# expandid about pole of order m 
# if c=infinity, then x is replaced by 1/y and expansion is around 0 
# for infinity, n=0,1,2,... 
#-------------------------------------------------------------------- 
local laurent_coeff::static:=proc(_self, 
        f, 
        x::symbol, 
        c, 
        m::integer, 
        n::integer,$) 
 
    local the_coeff,y::nothing,fy; 
 
    if c = infinity then 
        fy := eval(f,x=1/y); 
        if m-n>0 then 
           the_coeff := limit( diff( y^m * fy,y$(m-n)),y=0,right)/(m-n)!; 
        elif m-n=0 then 
            the_coeff := limit( y^m * fy,y=0,right)/(m-n)!; 
        else 
           the_coeff := 0; 
        fi; 
    else 
        if m-n>0 then 
            the_coeff := limit( diff( (x-c)^m * f,x$(m-n)),x=c,right)/(m-n)! ; 
        elif m-n=0 then 
           the_coeff := limit( (x-c)^m * f,x=c,right)/(m-n)!; 
        else 
           the_coeff := 0; 
        fi; 
    fi; 
 
    the_coeff := eval(the_coeff,[csgn=1,signum=1]); 
 
    return the_coeff; 
 
end proc; 
 
end module;

References

[1] Jerald J. Kovacic. An Algorithm for Solving Second Order Linear Homogeneous Diﬀerential Equations. J. Symb. Comput., 2(1):3–43, 1986.

[2] Carolyn J. Smith. A DISCUSSION AND IMPLEMENTATION OF KOVACICS ALGORITHM FOR ORDINARY DIFFERENTIAL EQUATIONS. Research Report CS-84-35, October 1984.

[3] B. David Saunders. An Implementation of Kovacic’s Algorithm for Solving Second Order Linear Homogeneous Diﬀerential Equations. Proceedings of the 1981 ACM Symposium on Symbolic and Algebraic Computation, 1981.

¹The complete Kovacic package in one mpl ﬁle accompany the arXiv version of this paper

²Wikipedia deﬁnes Liouvillian function as function of one variable which is the composition of a ﬁnite number of arithmetic operations $(+, -, \times, \div)$ , exponentials, constants, solutions of algebraic equations (a generalization of nth roots), and antiderivatives. Kovacic in his original paper says “Such a solution may involve exponentials, indeﬁnite integrals and solutions of polynomial equations. (As we are considering functions of a complex variable, we need not explicitly mention trigonometric functions, they can be written in terms of exponentials. Note that logarithms are indeﬁnite integrals and hence are allowed.”