This description is based on KOVACIC 1985 paper and not based on the Saunders paper.
We are given \(y^{\prime \prime }=ry\). It is assumed that the necessary conditions for case 1 have been met as given in the table above and \(r=\frac {s}{t}\) where \(\gcd \left ( s,t\right ) =1\) (in Maple this is done using the normal() command). The first step is to find the poles of \(r\) and the order of each pole. If there are no poles, then let the set of poles \(\Gamma \) will be empty.
If a pole \(x=c\) is of order 1 which means there is a factor \(\frac {1}{\left ( x-c\right ) }\) in the partial fractions decomposition of \(r\), then let \begin {align*} \left [ \sqrt {r}\right ] _{c} & =0\\ \alpha _{c}^{+} & =1\\ \alpha _{c}^{-} & =1 \end {align*}
If the pole \(c\) is of order 2, which means there is a factor \(\frac {1}{\left ( x-c\right ) ^{2}}\) in the partial fractions decomposition of \(r\), then let \begin {align*} \left [ \sqrt {r}\right ] _{c} & =0\\ \alpha _{c}^{+} & =\frac {1}{2}+\frac {1}{2}\sqrt {1+4b}\\ \alpha _{c}^{-} & =\frac {1}{2}-\frac {1}{2}\sqrt {1+4b} \end {align*}
Where \(b\) is the coefficient of \(\frac {1}{\left ( x-c\right ) ^{2}}\) in the partial fraction decomposition of \(r\). For example, if \(r=\frac {3}{\left ( x-2\right ) ^{2}}\), then \(x=2\) is a pole of order \(2\) and \(b=3\). The coefficients are found using undetermined coefficients method. (Examples below show how).
If the pole is of order \(4\) or \(6\) or \(8\) and so on, then it is a little bit more complicated. We write \(2v=order\). For example, if the pole was order 4, then \(v=2\) and if the pole was order 6, then \(v=3\) and so on. Notice that for case 1, which we are discussing here, if pole is of order larger than 2, then only poles of order \(4,6,8,\cdots \) are allowed. This is from the necessary condition. In this case, we add all terms involving \(\frac {1}{\left ( x-c\right ) ^{i}}\) for \(2\leq i\leq v\) in the Laurent series expansion of \(\sqrt {r}\) (not \(r\)) as follows\begin {align} \left [ \sqrt {r}\right ] _{c} & =\sum _{i=2}^{v}\frac {a_{i}}{\left ( x-c\right ) ^{i}}\tag {1}\\ & =\frac {a_{2}}{\left ( x-c\right ) ^{2}}+\frac {a_{3}}{\left ( x-c\right ) ^{3}}+\cdots +\frac {a_{v}}{\left ( x-c\right ) ^{v}}\nonumber \end {align}
For an example if the pole was of order \(6\), then \(v=3\). Therefore we need to add all terms in the Laurent series expansion of \(\sqrt {r}\) from \(v=3\) down to \(2\). As follows\begin {align*} \left [ \sqrt {r}\right ] _{c} & =\sum _{i=2}^{3}\frac {a_{i}}{\left ( x-c\right ) ^{i}}\\ & =\frac {a_{2}}{\left ( x-c\right ) ^{2}}+\frac {a_{3}}{\left ( x-c\right ) ^{3}} \end {align*}
Lets look at an example of the above before going to the next step. Assume \[ r=\frac {4x^{6}-8x^{5}+12x^{4}+4x^{3}+7x^{2}-20x+4}{4x^{4}}\] There is only one pole at \(x=0\) of order \(4\). Hence \(v=2\). We need to find the Laurent series of \(\sqrt {r}\) expanding around the specific pole \(c\) of order \(2v\). In Maple this is done using \(\operatorname {series}(\sqrt {r},x=c)\). \begin {equation} \sqrt {r}\approx \frac {1}{x^{2}}-\frac {5}{2}\frac {1}{x}-\frac {9}{4}-\frac {41}{8}x-\frac {443}{32}x^{2}+\cdots \tag {2} \end {equation} Hence \begin {align} \left [ \sqrt {r}\right ] _{c} & =\sum _{i=2}^{2}\frac {a_{i}}{\left ( x-c\right ) ^{i}}\nonumber \\ & =\frac {a_{2}}{\left ( x-c\right ) ^{2}}\nonumber \\ & =\frac {a_{2}}{\left ( x-0\right ) ^{2}} \tag {3} \end {align}
Comparing the above to Eq. 2 shows that the coefficient is \(a_{2}\) is (written now as just \(a\) to make it match the paper and use it in the following equation later on) \[ a=1 \] So the term \(a\) is the coefficient of \(\frac {a_{v}}{\left ( x-c\right ) ^{v}}\) in the Laurent series expansion of \(\sqrt {r}\) around \(x=c\). In implementation of the algorithm the method of undetermined coefficients is used instead of actually finding Laurent series for \(\sqrt {r}\) at \(x=c\).
Now that we found \(\left [ \sqrt {r}\right ] _{c}\) for poles \(>2\), we need to find its \(\alpha _{c}^{+},\alpha _{c}^{-}\) also. In this case \(\alpha _{c}^{+}=\frac {1}{2}\left ( \frac {b}{a}+v\right ) \) and \(\alpha _{c}^{-}=\frac {1}{2}\left ( -\frac {b}{a}+v\right ) \). Where \(a\) is the one we just found above. But what is \(b\) here? \(b\) is the coefficient of the term \(\frac {1}{\left ( x-c\right ) ^{v+1}}\) in \(r\) minus the coefficient of \(\frac {1}{\left ( x-c\right ) ^{v+1}}\) in \(\left [ \sqrt {r}\right ] _{c}\) which we found above in (2). For an example, using the above \(r\) its Laurent series expansion around \(x=0\) is\[ r\approx x^{2}-2x+3+\frac {1}{x}+\frac {7}{x^{2}}-\frac {5}{x^{3}}+\frac {1}{x^{4}}\] Then since \(a=1\) from earlier and since \(v=2\) here (since pole of order 4) then we look above for the coefficient of the term \(\frac {1}{\left ( x-0\right ) ^{v+1}}=\frac {1}{\left ( x-0\right ) ^{3}}\) in \(r\) itself. We see this is \(-5\). Now we need to subtract from this value the coefficient of \(\frac {1}{\left ( x-0\right ) ^{v+1}}=\frac {1}{\left ( x-0\right ) ^{3}}\) from \(\left [ \sqrt {r}\right ] _{c}\) series from Eq (3). But since \(\left [ \sqrt {r}\right ] _{c}=\frac {1}{\left ( x-0\right ) ^{2}}\) then there is no term \(\frac {1}{\left ( x-0\right ) ^{3}}\). Which means \begin {align*} b & =-5-0\\ & =-5 \end {align*}
Therefore for this example\begin {align*} \alpha _{c}^{+} & =\frac {1}{2}\left ( \frac {b}{a}+v\right ) =\frac {1}{2}\left ( -\frac {5}{1}+2\right ) =-\frac {3}{2}\\ \alpha _{c}^{-} & =\frac {1}{2}\left ( -\frac {b}{a}+v\right ) =\frac {1}{2}\left ( \frac {5}{1}+2\right ) =\frac {7}{2} \end {align*}
We are now done with finding everything we need related to poles. The above needs to be done for each pole \(c\) in \(r\).
We see that for each pole, we need to calculate 3 items. They are \(\left [ \sqrt {r}\right ] _{c},\alpha _{c}^{+},\alpha _{c}^{-}\).
Now we switch attention to the \(O\left ( \infty \right ) \) order. This is much easier. This is the order of \(r=\frac {s}{t}\) at infinity which is found from \(\deg \left ( t\right ) -\deg \left ( s\right ) \). There are also three cases to consider.
If \(O\left ( \infty \right ) >2\) then we write \begin {align*} \left [ \sqrt {r}\right ] _{\infty } & =0\\ \alpha _{\infty }^{+} & =0\\ \alpha _{\infty }^{-} & =1 \end {align*}
If \(O\left ( \infty \right ) =2\) then \(\left [ \sqrt {r}\right ] _{\infty }=0\). Now we calculate \(b\) for this case. This is given by the leading coefficient of \(s\) divided by the leading coefficient of \(t\) when \(\gcd \left ( s,t\right ) =1\). In this case \begin {align*} \left [ \sqrt {r}\right ] _{\infty } & =0\\ \alpha _{\infty }^{+} & =\frac {1}{2}+\frac {1}{2}\sqrt {1+4b}\\ \alpha _{\infty }^{-} & =\frac {1}{2}-\frac {1}{2}\sqrt {1+4b} \end {align*}
Where here \(b\) is the coefficient of \(\frac {1}{x^{2}}\) in the Laurent series expansion of \(r\) at \(\infty \). But we do not need to find Laurent series expansion of \(r\) at \(\infty \) to find \(b\) here. It can be found using \(b=\frac {lcoeff\left ( s\right ) }{lcoeff\left ( t\right ) }\) where \(r=\frac {s}{t}\) and \(\gcd \left ( s,t\right ) =1\). And \(lcoeff\) is the leading coefficient. For example, if \(r=\frac {1+5x}{2x^{2}}\) then \(b=\frac {1}{2}\). If we took the Laurent series of \(r\) at \(\infty \) which in Maple can be done using the command \(series(r,x=\infty )\) then we will get \(\frac {5}{2x}+\frac {1}{2}\frac {1}{x^{2}}\) which also give \(b=\frac {1}{2}\).
And finally, if \(O\left ( \infty \right ) \) \(\leq 0\), then \(O\left ( \infty \right ) \) has to be negative and even number (conditions for case 1). Let the order of \(r\) at \(\infty \) be \(-2v\leq 0\). Then now \(\left [ \sqrt {r}\right ] _{\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 \).\[ \left [ \sqrt {r}\right ] _{\infty }=ax^{v}+z_{1}x^{v-1}+\cdots +z_{n}\] And \(b\) is the coefficient of \(x^{v-1}\) in \(r\) minus the coefficient of \(x^{v-1}\) in \(\left ( \left [ \sqrt {r}\right ] _{\infty }\right ) ^{2}\). Then \begin {align*} \alpha _{\infty }^{+} & =\frac {1}{2}\left ( \frac {b}{a}-v\right ) \\ \alpha _{\infty }^{-} & =\frac {1}{2}\left ( -\frac {b}{a}-v\right ) \end {align*}
This completes step 1 of the algorithm. We have found \(\left [ \sqrt {r}\right ] _{c}\) for each pole and associated \(\alpha _{c}^{+},\alpha _{c}^{-}\) and also found \(\left [ \sqrt {r}\right ] _{\infty }\) and its associated \(\alpha _{\infty }^{+},\alpha _{\infty }^{-}\). So, what will we do with these? In step 2 these are used to find all possible values of what is called \(d\). For each non negative \(d\), we will find a candidate \(\omega \). And use this candidate \(\omega \) to find \(P\left ( x\right ) \) by solving \(P^{\prime \prime }+2\omega P^{\prime }+\left ( \omega ^{\prime }+\omega ^{2}-r\right ) P=0\) (linear algebra problem). If we are able to find a \(P\left ( x\right ) \) for any one candidate \(\omega \) then we stop and we have found the solution \(y=p\left ( x\right ) e^{\int \omega dx}\) to the \(y^{\prime \prime }=ry\). Examples below will show how all this works.
Recall that from step 1 we have found \(\left [ \sqrt {r}\right ] _{c}\) and its associated \(\alpha _{c}^{+},\alpha _{c}^{-}\) (this is done for each pole of \(r\)) and we have found \(\left [ \sqrt {r}\right ] _{\infty }\) and its associated \(\alpha _{\infty }^{+},\alpha _{\infty }^{-}\). From these we now found a possible \(d\) values and trying each \(d\geq 0\). The value of \(d\) is found using the following for each combination of \(s\left ( c\right ) \) where \(s\left ( c\right ) \) is \(+\) or \(-\) \[ d=\alpha _{\infty }^{\pm }-\sum _{c}\alpha _{c}^{\pm }\] We keep only the non negative values of \(d\). It is important to note that we have to find an integer positive value for \(d\) to continue. If no such value is found from the above, then we stop here as this means no Liouvillian solution exist using case 1. Then we go to case two or case three if it is available.
If we do find \(d\geq 0\), then we now find corresponding candidate \(\omega _{d}\) using\begin {align*} \omega _{d}&=\sum _{c}\left ( (\pm ) \left [ \sqrt {r}\right ] _{c}+\frac {\alpha _{c}^{\pm }}{x-c}\right ) + (\pm ) \left [\sqrt {r}\right ]_{\infty } \end {align*}
In this step we first attempt to find a polynomial \(p\left ( x\right ) \) of degree \(d\), for \(\omega \) found in step 2. This is done by solving \[ p^{\prime \prime }+2\omega p^{\prime }+\left ( \omega ^{\prime }+\omega ^{2}-r\right ) p=0 \] For example, if \(d=2\), then we let \(p\left ( x\right ) =x^{2}+ax+b\) and if \(\omega \) happened to be say \(\frac {1}{x^{2}}-\frac {3}{2x}+x-1\), then by substituting these in the above, we can solve for \(a,b\) (if a solution exist). Then the solution to \(y^{\prime \prime }=ry\) is \(y=p\left ( x\right ) e^{\int \omega dx}\). If the degree \(d=1\) then we guess \(p\left ( x\right ) =x+a\) and try to solve for \(a\). If the degree \(d=0\), then we let \(p\left ( x\right ) =1\), a constant. In the special case of \(p\left ( x\right ) =1\), there is no coefficients \(a_{i}\) to solve for. So we would just need to verify that \[ \omega ^{\prime }+\omega ^{2}-r=0 \] In this case.
This completes the full algorithm for case 1. We will now go over many examples for case 1, showing how to implement this algorithm for each example.
The hardest part of the kovacic algorithm is just finding all the \(\left [ \sqrt {r}\right ] _{c},\alpha _{c}^{\pm },\left [ \sqrt {r}\right ] _{\infty },\alpha _{\infty }^{\pm }\). Once these are found, the rest of the algorithm is much more direct.
The following diagram summarized the above for case one.