On the Solutions of Second-Order Differential Equations with Polynomial Coefﬁcients: Theory, Algorithm, Application

: The analysis of many physical phenomena is reduced to the study of linear differential equations with polynomial coefﬁcients. The present work establishes the necessary and sufﬁcient conditions for the existence of polynomial solutions to linear differential equations with polynomial coefﬁcients of degree n , n − 1, and n − 2 respectively. We show that for n ≥ 3 the necessary condition is not enough to ensure the existence of the polynomial solutions. Applying Scheffé’s criteria to this differential equation we have extracted n generic equations that are analytically solvable by two-term recurrence formulas. We give the closed-form solutions of these generic equations in terms of the generalized hypergeometric functions. For arbitrary n , three elementary theorems and one algorithm were developed to construct the polynomial solutions explicitly along with the necessary and sufﬁcient conditions. We demonstrate the validity of the algorithm by constructing the polynomial solutions for the case of n = 4. We also demonstrate the simplicity and applicability of our constructive approach through applications to several important equations in theoretical physics such as Heun and Dirac equations.


Introduction
Differential equations with polynomial coefficients have played an important role not only in understanding engineering and physics problems , but also as a source of inspiration for some of the most crucial results in special functions and orthogonal polynomials [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. The classical differential equation α 2 r 2 + α 1 r + α 0 y (r) + β 1 r + β 0 y (r) − τ 0 y(r) = 0 , (1) has been a cornerstone of many fundamental results in special functions since the early work of L. Euler [33], see also References [47][48][49]. More recently, the Heun-type equation [41][42][43], r α 3 r 2 + α 2 r + α 1 y (r) + β 2 r 2 + β 1 r + β 0 y (r) − τ 1 r + τ 0 y(r) = 0 , has become a classic equation in mathematical physics. These two equations are members of the general differential equation P n (r) y (r) + Q n−1 (r) y (r) − R n−2 (r) y(r) = 0 , n ≥ 2 , P n (r) = n ∑ k=0 α k r k , Q n−1 (r) = n−1 ∑ k=0 β k r k , R n−2 (r) = where {α k } n k=0 , {β k } n−1 k=0 , {τ k } n−2 k=0 are real parameters. We assume that the leading term of R n−2 = 0 and at least one of the leading terms α n and β n−1 is not zero. The purpose of the present work is to provide an answer to the following question: "Under what conditions on the equation parameters {α k } n k=0 , {β k } n−1 k=0 , and {τ k } n−2 k=0 does the differential Equation With the general theory of linear differential equations as a guide [50,51], the zeros of the leading polynomial P n (r) classify the solutions of the Equation (3). The point r = ζ is called an ordinary point if Q n−1 (r)/P n (r) and R n−2 (r)/P n (r) are analytic functions at r = ζ, or a regular singular point if Q n−1 (r)/P n (r) and R n−2 (r)/P n (r) are not analytic at r = ζ but the products (r − ζ) Q n−1 (r)/P n (r) and (r − ζ) 2 R n−2 (r)/P n (r) are analytic at that point. Using this characterization of singularities, S. Bochner [48] classified the polynomial solutions of (1) in terms of the classical orthogonal polynomials. A different approach which depends on the parameters of the leading coefficient P n (r), was introduced in Reference [13] to study the polynomial solutions of (1). It was shown that there are seven possible nonzero leading polynomials depending on the combination of the nonzero α i parameters, for i = 0, 1, 2. By analyzing each of these cases, the authors [13] were able to explicitly construct all the polynomial solutions of Equation (1) in terms of hypergeometric functions along with the associated weight functions.
The nth-degree leading polynomial coefficient P n (r) of the differential Equation (3) contains (2 n+1 − 1) nonzero polynomials that depend on the nonzero values of α i , i = 0, 1, . . . , n. Out of this polynomial set, there are 2 n differential equations for which r = 0 is an ordinary point, (2 n−1 + 2 n−2 ) differential equations for which r = 0 is a regular singular point, and the remaining (2 n−2 − 1) differential equations with irregular singular points that fall outside of the scope of this present work.
The present work establishes polynomial solutions in simple and constructive forms to (3) along with all necessary and sufficient conditions they require. This work is done independently from that of Theorem 1, although can be verified using the Asymptotic Iteration Method.
The paper is organized as follows: • In Section 2, we introduce the concept of a necessary but not sufficient condition to find polynomial solutions, and then demonstrate how the inverse square-root potential [22][23][24][25] is one such example that does not have polynomial solutions yet satisfies the necessary condition.

•
In Section 3, Scheffé criteria is revised to analyze Equation (3) and to generate exactly solvable equations where the coefficients of their power series solutions are easily computed using a two-term recurrence relation. The closed form solutions of these equations are written in terms of the generalized hypergeometric functions where the Pochhammer symbol (α) k is defined in terms of the Gamma function as In Section 4 we present three theorems corresponding to the cases: (1) α 0 = 0, (2) α 0 = 0, α 1 = 0, and (3) α 0 = α 1 = β 0 = 0, α 2 = 0 which are required to establish the polynomial solutions of Equation (3). We shall show that for the mth-degree polynomial solutions, there are n + m − 1 conditions that ultimately assemble these polynomials. At the end of the section we briefly discuss the Mathematica R program used to generate the solutions of these equations. • Lastly, Section 5 demonstrates the validity of these constructions through the application of our results to some problems that have appeared in mathematics and physics literature.

A Necessary but Not Sufficient Condition
We begin by stating the necessary condition for the existence of polynomial solutions to the general differential equation given by (3): Theorem 2. A necessary condition for the differential Equation (3) to have a polynomial solution of degree m is Proof. Substitute the polynomial solution y = ∑ m i=0 C i r i into the differential Equation (3). After collecting the common terms of the common r power, the necessary condition for the leading coefficient C m to be nonzero is τ n−2 = α n m (m − 1) + β n−1 m for all m ≥ 0 follows.
For n = 2, Equation (7) represents both the necessary and sufficient conditions for the polynomial solutions of (3). Although the condition (7) is necessary, it is not sufficient for n ≥ 3 to guarantee the polynomial solutions of (3). There are precisely n − 2 additional conditions that are sufficient to ensure the existence of such polynomials. These additional conditions can be understood as constraints that relate the remaining parameters of R n−2 (r) with the coefficients {α k } n k=0 and {β k } n−1 k=0 in P n (r) and Q n−1 (r), respectively. In later sections, we shall devise a procedure to find these sufficient conditions and provide a method of evaluating the corresponding polynomial solutions explicitly.
The remainder of the section investigates the inverse square-root potential. We shall show that despite satisfying the necessary condition in Theorem 2 that the Schrödinger equation with this potential as mentioned in the literature has no polynomial solutions, thus demonstrating that the condition (7) is necessary but not sufficient.

The Inverse Square-Root Potential
The exact solutions of the radial Schrödinger equation were recently studied by a number of authors [21][22][23][24][25]. The differential Equation (8) has two singular points, one regular at r = 0 and another irregular at r = ∞ of rank 3 [50]. To analyze these singular points further, we employ a change of variables z = √ 2 κ r, for some constant κ to be determined shortly. Equation (8) may then be expressed as As z → 0, the solution to the differential Equation (9) behaves asymptotically as the solution of the Euler equation, z 2 ψ (z) − z ψ (z) − 4 l (l + 1) ψ(z) = 0, which has a physically acceptable solution ψ(z) ∼ z 2 (l+1) . Meanwhile as z → ∞, the solution to the differential equation behaves asymptotically as the solution of To examine the solution of this equation we complete the square to The substitution x = ((v/ √ 2 κ 3 ) − z) allows us to compare the solution of the equation to the solvable Schröinger equation with the harmonic oscillator potential, as z → ∞, via Hence, we may assume the solution of Equation (9) to be up to a normalization constant N. Upon substituting the ansatz solution (12) into Equation (9), we obtain a second-order differential equation (an example of a biconfluent Heun equation [41]) for f (z) as subject to κ = √ −E. Clearly, Equation (13) is of the form of Equation (3) with n = 3 with the parameters Using Theorem 2, for the existence of an mth-degree polynomial solution it is necessary that Note that this result can be also deduced using the Bethe Ansatz method [12,54]. On the other hand, an application of the Frobenius method establishes the following three-term recurrence relation for the coefficients of the polynomial solutions f m (r) = ∑ m j=0 C j r j : where j = 0, 1, 2, . . . , m, m + 1 and C j = 0 for j < 0. Implementing the necessary condition (14), with λ = − √ 4 + 4l + 2m < 0, Equation (15) then reads The consistency of the linear homogeneous system generated by the three-term recurrence relation (16), for j = 0, 1, · · · , m + 1, enforces the vanishing of the determinant, denoted by ∆ (m+1) , that establish the sufficient condition Therefore, if the determinant ∆ m+1 is different than zero, the only possible solution to the linear system is the trivial solution (i.e., the zero solution). Hence any nonzero value of the determinant ∆ m+1 is an indication that a polynomial solution is not possible, which is the case for the coefficients of Equation (13).

Remark 1.
The point here is not solving the Schrödinger equation with the inverse-root potential; this is already done in several manuscripts [21][22][23][24][25]. We stress that the necessary condition is not enough to guarantee polynomial solutions. Also, the necessary and sufficient constraints must have common roots that are physically acceptable.

Remark 2.
In Reference [55], the author claims the existence of polynomial solutions of a differential equation which matches (13) that satisfy the necessary condition (7). The author missed the possibility of the non-vanishing determinant (4), in his work, for some β 0 and β 2 .

Remark 3.
For the applications in physics, there is nothing special about the necessity of the polynomial solutions. However, the problem of the existence of polynomial solutions is important. Studying the problem in its general form and developing mathematical tools to treat it is significant on its own.

Scheffé's Criteria: Two-Term Recurrence Relation
Generally speaking, recurrence relations with more than two terms are difficult to solve explicitly. Differential equations that are known to have a two-term recurrence relation for their power series solutions guarantee the solvability of such equations.
In a paper presented to the American Mathematical Society (1941), H. Scheffé [34] establishes criteria for the necessary and sufficient conditions for differential equations of the form to have a two-term recurrence relation. Here p j (r) = 0, j = 0, 1, 2, are analytic functions (not necessarily polynomials) in some region consisting of all points in an arbitrary neighbourhood of a regular singular point r = r 0 except the point r 0 itself. We adopt this criterion to examine the differential Equation (18) with and provide a formula for the two-term recurrence relation. Without loss of generality, we shall take the ordinary or singular point as r = 0, otherwise a simple shifting of r is applied first.
Theorem 3. The necessary and sufficient conditions for the differential Equation (18) to have a two-term recurrence relation between successive coefficients in its series solution, relative to the ordinary or regular singular point r = 0, is that in the neighbourhood of r = 0 Equation (18) can be written as: where, for m ∈ Z , h ∈ Z + , and at least one of the product q j,0 q i,h , for i, j = 0, 1, 2, is nonzero. The two-term recurrence formula is given by where λ = λ 1 , λ 2 are the roots of the indicial equation The closed form of the series solution can be written in terms of the generalized hypergeometric function as As an example, we consider the differential equation mentioned in the introduction. We want to extract the possible equations that have two-term recurrence relations between successive coefficients in their series solutions. Comparing Equation (25) with (20), we find that h + 2 − m = 2 and h + 2 − m = 2 > 2 − m ≥ 0, which leads to h = m, 2 > 2 − m ≥ 0. Consequently 0 < m ≤ 2 which gives us two cases: (m, h) = (1, 1), and (m, h) = (2, 2).

•
In the case of m = 1 and h = 1, Equation (20) reads: and the corresponding Equation (25) with q 2,0 = α 1 , q 2,h=1 = α 2 , q 1,0 = β 0 , q 1,h=1 = β 1 , q 0,h=1 = γ 0 and q 0,0 = 0 reads The two-term recurrence formula for the solution of Equation (27) is where λ are the zeros of the indicial equation The closed form solutions in terms of the Gauss hypergeometric functions, respectively, are: and • In this case of m = 2 and h = 2, Equation (20) reads and the corresponding Equation (25) reads The two-term recurrence relation for the solutions of Equation (32) reads where, λ are the roots of the indicial equation λ (λ − 1) = 0, α 0 = 0, with closed form solutions and u 2 (r, λ = 1) Remark 4. For the admissible values of the equation parameters, the recurrence relations (28) and (33) are generic formulas for the two-term recurrence relations for the series solutions of the differential Equations (27) and (32) respectively. That is to say, by assigning admissible values of the parameters, (27) and (32) will generate solvable differential equations.

Remark 5.
Polynomial solutions for the differential Equations (27) and (32) can be easily obtained using suitable values of the equation parameters to terminate the recurrence relations (28) and (33).

Theorems and Algorithm
In this section we present three elementary theorems that classify the polynomial solutions of the general differential Equation (3). Following that, we will give a brief description of the Mathematica R program that was developed to accompany the present work. A link to the complete program, available for direct use, is also provided for direct applications of these theorems.

Theorems
Theorem 4 applies to 2 n equations of (3) and gives a recurrence relation needed to compute polynomial solutions in the neighbourhood of an ordinary point. Theorems 5 and 6 give the polynomial solutions for the 2 n−1 + 2 n−2 equations from (3) about a regular singular point. The proofs of Theorems 4 and 5 may be found in the Appendix A, while the proof of Theorem 6 follows closely to the proof of Theorem 5 and is therefore omitted. For α 0 = 0, the second-order linear differential Equation (3) admits the solution y(r) = P m (r) = ∑ m i=0 C i r i in the neighbourhood of the ordinary point r = 0 and valid to the nearest nonzero singular point of P m (r) = 0. Here P m (r) is an mth-degree polynomial if for all = 0, 1, 2, . . . , m + n − 2, and τ i = β i = 0 if i < 0. The equation corresponding to the value of = m + n − 2 gives the necessary condition (7) for the existence of the mth-degree polynomial solution The remaining (m + n − 2) linear equations consist of (n − 2) sufficient conditions that relate {α k } n For α 0 = 0 and α 1 = 0, the second-order linear differential Equation (3) admits the solution y(r) = r s P m (r) = r s ∑ m i=0 C i r i in the neighbourhood of the regular singular point r = 0, where P m (r) is an mth-degree polynomial if n−1 ∑ j=0 α n−j ( − n + j + s + 2) ( − n + j + s + 1) + β n−j−1 ( − n + j + s + 2) − τ n−j−2 C −n+j+2 = 0 , (53) for all = 0, 1, 2, . . . , m + n − 2, and τ i = 0 if i < 0. Here s is a root of the indicial equation: s(s − 1) + (β 0 /α 1 )s = 0, so s ∈ 0 , 1 − (β 0 /α 1 ) .
For α 0 = α 1 = β 0 = 0 and α 2 = 0, the second-order linear differential Equation (3) admits the solution y(r) = r s P m (r) = r s ∑ m i=0 C i r i in the neighbourhood of the regular singular point r = 0, where P m (r) is an mth-degree polynomial if (54) for all = 0, 1, 2, . . . , m + n − 2. Here s is a nonzero root of the indicial equation: s(s − 1) + (β 0 /α 1 )s − for which the recurrence relation (54) reduces to: Clearly when = m and if C m = 0 then Equation (56) gives the necessary condition: as expected for the solutions of the Cauchy-Euler equation.

The Mathematica R Program
A Mathematica R algorithm and accompanying GUI ( Figure 1) were developed to assist with this work. The algorithm uses Cramer's rule to compute the coefficients of the polynomial solutions {C i } m i=1 for the Equation (3) in the case where the equation is in the form of Theorems 4-6. In addition to computing the coefficients, the program also generates the sufficient and necessary conditions. A complete Mathematica R program may be obtained by contacting the authors, or via the following Github repository: https://github.com/KyleBryenton/OPSOLDE.  [2,9], the coefficient field will resize appropriately and permit the user to choose the values for each coefficient. A working example is provided: by selecting "Dirac Coefficients" and the program will generate the results discussed in Section 5.2 of this article.

Applications
The simplicity of the above theorems lay in their constructive approach to explicitly computing the polynomial solution coefficients and to provide all the necessary and sufficient conditions for the existence of these polynomials. Further, they can be easily implemented using any available symbolic algebra system. In this section, we show how to make use of these theorems.

Conclusions
Given a linear differential equation with polynomial coefficients of order n, n − 1, and n − 2 respectively, this work provides a constructive and straightforward approach for finding all possible polynomial solutions along with the necessary and sufficient conditions. A Mathematica program which solves these differential equations for an integer value of n ∈ [2,9] is available for download (see Section 4.2), but the algorithms are also provided so that the reader may adopt them in their preferred computer software. Although there may be other sophisticated approaches available in the literature for solving differential equations [46], the method presented here has its simplicity characterized by recurrence relations. Aside from the examples discussed in this present work, the authors explored a large number of linear differential equations which appeared in mathematics and physics literature, and have found exact consistency with the results obtained by other researchers [3][4][5][6][7][8][9][10]43,55].