Quasi-Exactly Solvable Differential Models: A Canonical Polynomials Approach

Abstract

A quasi-exactly solvable model refers to any second order differential equation with polynomial coefficients of the form A(x)y’’(x)+B(x)y’(x)+C(x)y(x)=0 where a pair of exact polynomials {y(x), C(x)} with respective degrees {deg[y]=n, deg[C]=p} are to be found simultaneously in terms of the coefficients of two given polynomials {A(x), B(x)}. The existing methods for solving quasi-exactly solvable models require the solution of a system of nonlinear algebraic equations of which the dimensions depend on n, the degree of the exact polynomial solution y(x). In this paper, a new method employing a set of polynomials, called canonical polynomials, is proposed. This method requires solving a system of nonlinear algebraic equations of which the dimensions depend only on p, the degree of C(x), and do not vary with n. Several examples are implemented to testify the efficiency of the proposed method.

Share and Cite:

El-Daou, M. , AlZanki, T. and Al-Mutawa, N. (2019) Quasi-Exactly Solvable Differential Models: A Canonical Polynomials Approach. American Journal of Computational Mathematics, 9, 48-60. doi: 10.4236/ajcm.2019.92004.

1. Introduction

Let D be a linear differential operator of degree $\nu \ge 1$ with polynomial coefficients, and let $y\left(x\right)$ be an exact solution of the differential equation:

$\left(Dy\right)\left(x\right)\equiv \underset{i=0}{\overset{\nu }{\sum }}{A}_{i}\left(x\right)\frac{{\text{d}}^{i}y}{\text{d}{x}^{i}}={A}_{\nu +1}\left(x\right)$ (1)

where $\left\{{A}_{i}\left(x\right);i=0,1,\cdots ,\nu +1\right\}$ are polynomials in x for $i=0,1,\cdots ,\nu +1$.

Among the very efficient techniques to solve equations of the form (1) is the Tau Method. This method was invented by Lanczos in [1] to solve simple equations and it was extended later by Ortiz [2] to treat more sophisticated problems with applications in a wide range of scientific fields. The basic idea of the Tau Method consists of finding polynomial solutions to Equation (1) in terms of a special set of polynomials associated with D called canonical polynomials, denoted by $\left\{{Q}_{k}\left(x\right);k\ge 0\right\}$. These are defined as the pseudo inverse image of the set of the power functions $\left\{{x}^{k};k\ge 0\right\}$ by D.

In this paper, we use the canonical polynomials to develop a new method to solve quasi-exactly solvable (QES) second order ordinary differential equations. The QES problem, considered in this paper, can be stated as follows: Let p and n be two nonnegative integers, and suppose that $\left\{A\left(x\right),B\left(x\right),C\left(x\right)\right\}$ are polynomials with expressions

$A\left(x\right)=\underset{i=0}{\overset{p+2}{\sum }}{A}_{i}{x}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}B\left(x\right)=\underset{i=0}{\overset{p+1}{\sum }}{B}_{i}{x}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{and}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}C\left(x\right)=\underset{i=0}{\overset{p}{\sum }}{C}_{i}{x}^{i}.$

The coefficients $\left\{{C}_{i};i=0,1,\cdots ,p\right\}$ are found in terms of $\left\{{A}_{i};i=0,1,\cdots ,p+2\right\}$ and $\left\{{B}_{i};i=0,1,\cdots ,p+1\right\}$ so that the function $y\left(x\right)$ defined implicitly by the differential equation

$\left(Dy\right)\left(x\right)=A\left(x\right)\frac{{\text{d}}^{2}y}{\text{d}{x}^{2}}+B\left(x\right)\frac{\text{d}y}{\text{d}x}+C\left(x\right)y\left(x\right)=0$ (2)

is an exact polynomial solution of degree n and written as $y\left(x\right)={\sum }_{i=0}^{n}\text{ }{y}_{i}{x}^{i}$, where $\left\{{y}_{i};i=0,1,2,\cdots ,n\right\}$ are determined in terms of $\left\{{A}_{i},{B}_{i},{C}_{i};i\ge 0\right\}$.

Quasi-exactly solvable problems have applications in engineering, chemistry and physics among many other fields. This includes a wide range of mathematical settings that involve Schrödinger equations describing problems in quantum mechanics, for example, anharmonic singular potentials, coulombically repelling electrons on a multidimensional sphere, Planer Dirac electron in magnetic fields and Kink stability analysis among many other problems (see [3] [4] [5] [6] and the references given therein).

QES problems are usually solved by means of the Functional Bethe Ansatz method or by a constraint polynomial approach (see [7] ). With the latter, we derive a $n×n$ system of algebraic equations, which gives the n roots of $y\left(x\right)$, and afterward the p coefficients of $C\left(x\right)$ are calculated. If one wishes to increase the order n of polynomial $y\left(x\right)$, then a new nonlinear algebraic system of higher dimensions has to be resolved. With our proposed method that employs the canonical polynomials, we solve first a $p×p$ system of algebraic equations giving the coefficients of $C\left(x\right)$ and then the coefficients of $y\left(x\right)$ are computed. With this method, a polynomial solution $y\left(x\right)$ of higher order n can be obtained without increasing dimensions $p×p$ of the algebraic system.

Section 2 introduces the canonical polynomials in the context that fits the problem under consideration. In Section 3 the main results of this paper are given and discussed. Sections 4 and 5 illustrate the implementation of the proposed method.

2. The Canonical Polynomials

In this section we recall the main features of the canonical polynomials associated with D (see Ortiz [2] ). We begin with this definition.

Definition 1. For any integer $k\ge 0$,

1) ${P}_{k}\left(x\right)$ is called the kth generating polynomial of D if ${P}_{k}\left(x\right)=D{x}^{k}$.

2) ${Q}_{k}^{*}\left(x\right)$ is called a kth canonical function of D if ${Q}_{k}^{*}\left(x\right)$ satisfies the differential equation $D{Q}_{k}^{*}={x}^{k}$.

In [2] , Ortiz presented an algorithm in the form of a self-starting recursive formula for computing the ${Q}_{k}^{*}$ 's associated with linear differential operators of arbitrary degree $\nu \ge 1$. The following theorem gives a variation of this algorithm for second order differential operators of the form (2) in which this paper is concerned:

Theorem 1. The canonical functions associated with the differential operator (2) can be generated by the following recursion:

${Q}_{k+p}^{*}=\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(k\right)}{Q}_{k+i-2}^{*}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}k\ge 0$ (3)

where ${f}_{i-2}^{\left(k\right)}:=k\left(k-1\right){A}_{i}+k{B}_{i-1}+{C}_{i-2}$ for any $i\ge 0$ with the convention that ${A}_{i}={B}_{i}={C}_{i}=0$ for $i<0$.

Proof. For any $k\ge 0$, let us expand the kth generating polynomial ${P}_{k}\left(x\right)$ in terms of $\left\{{x}^{i};i=0,1,2,\cdots \right\}$ :

$\begin{array}{l}{P}_{k}\left(x\right):=D{x}^{k}=A\left(x\right)k\left(k-1\right){x}^{k-2}+B\left(x\right)k{x}^{k-1}+C\left(x\right){x}^{k}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }=k\left(k-1\right)\underset{i=0}{\overset{p+2}{\sum }}{A}_{i}{x}^{k+i-2}+k\underset{i=0}{\overset{p+1}{\sum }}{B}_{i}{x}^{k+i-1}+\underset{i=0}{\overset{p}{\sum }}{C}_{i}{x}^{k+i}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }=k\left(k-1\right){A}_{p+2}{x}^{k+p}+k{B}_{p+1}{x}^{k+p}+{C}_{p}{x}^{k+p}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }\text{\hspace{0.17em}}\text{\hspace{0.17em}}+k\left(k-1\right)\underset{i=0}{\overset{p+1}{\sum }}{A}_{i}{x}^{k+i-2}+k\underset{i=0}{\overset{p}{\sum }}{B}_{i}{x}^{k+i-1}+\underset{i=0}{\overset{p-1}{\sum }}{C}_{i}{x}^{k+i}\end{array}$

$\begin{array}{l}=\left[k\left(k-1\right){A}_{p+2}+k{B}_{p+1}+{C}_{p}\right]{x}^{k+p}+k\left(k-1\right)\underset{i=0}{\overset{p+1}{\sum }}{A}_{i}{x}^{k+i-2}\\ \text{ }\text{\hspace{0.17em}}\text{ }\text{ }+k\underset{i=0}{\overset{p+1}{\sum }}{B}_{i-1}{x}^{k+i-2}+\underset{i=0}{\overset{p+1}{\sum }}{C}_{i-2}{x}^{k+i-2}\\ =\left[k\left(k-1\right){A}_{p+2}+k{B}_{p+1}+{C}_{p}\right]{x}^{k+p}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }+\underset{i=0}{\overset{p+1}{\sum }}\left[k\left(k-1\right){A}_{i}+k{B}_{i-1}+{C}_{i-2}\right]{x}^{k+i-2}\\ \equiv {f}_{p}^{\left(k\right)}{x}^{k+p}+\underset{i=0}{\overset{p+1}{\sum }}{f}_{i-2}^{\left(k\right)}{x}^{k+i-2}.\end{array}$ (4)

Since, by Definition 1, $D{Q}_{k+i-2}^{*}={x}^{k+i-2}$, and since D is a linear operator, (4) becomes:

$D{x}^{k}={f}_{p}^{\left(k\right)}{x}^{k+p}+D\left\{\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(k\right)}{Q}_{k+i-2}^{*}\right\},$ (5)

which implies that

$D\left\{{x}^{k}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(k\right)}{Q}_{k+i-2}^{*}\right\}={f}_{p}^{\left(k\right)}{x}^{k+p},$

from which follows the desired recursion (3):

${Q}_{k+p}^{*}=\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(k\right)}{Q}_{k+i-2}^{*}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\right)$ $\square$

For instance, when $k=0$, Equation (3) gives

$\begin{array}{c}{Q}_{p}^{*}=\frac{1}{{C}_{p}}\left\{1-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{C}_{i-2}{Q}_{i-2}^{*}\right\}=\frac{1}{{C}_{p}}\left\{1-\underset{i=0}{\overset{p-1}{\sum }}\text{ }{C}_{i}{Q}_{i}^{*}\right\}\\ =\frac{1}{{C}_{p}}-\frac{1}{{C}_{p}}\left({C}_{0}{Q}_{0}^{*}+{C}_{1}{Q}_{1}^{*}+\cdots +{C}_{p-1}{Q}_{p-1}^{*}\right)\equiv {Q}_{p}\left(x\right)+{R}_{p},\end{array}$

where

${Q}_{p}\left(x\right):=\frac{1}{{C}_{p}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{\hspace{0.17em}}{R}_{p}:=-\frac{1}{{C}_{p}}\left({C}_{0}{Q}_{0}^{*}+{C}_{1}{Q}_{1}^{*}+\cdots +{C}_{p-1}{Q}_{p-1}^{*}\right).$ (6)

That is, ${Q}_{p}^{*}$ has a polynomial component ${Q}_{p}\left(x\right)$ of degree 0 and a residual component ${R}_{p}\in \text{span}\left\{{Q}_{i}^{*};i=0,1,\cdots ,p-1\right\}$, which is the subspace generated by $\left\{{Q}_{0}^{*},{Q}_{1}^{*},\cdots ,{Q}_{p-1}^{*}\right\}$.

For $k=1$, we have

$\begin{array}{c}{Q}_{p+1}^{*}=\frac{1}{{B}_{p+1}+{C}_{p}}\left\{x-\underset{i=0}{\overset{p+1}{\sum }}\left({B}_{i-1}+{C}_{i-2}\right){Q}_{i-1}^{*}\right\}\\ =\frac{1}{{B}_{p+1}+{C}_{p}}\left\{x-\left({B}_{p}+{C}_{p-1}\right){Q}_{p}^{*}-\underset{i=0}{\overset{p}{\sum }}\left({B}_{i-1}+{C}_{i-2}\right){Q}_{i-1}^{*}\right\}\\ =\frac{1}{{B}_{p+1}+{C}_{p}}\left\{x-\left({B}_{p}+{C}_{p-1}\right)\left(\frac{1}{{C}_{p}}+{R}_{p}\right)-\underset{i=0}{\overset{p}{\sum }}\left({B}_{i-1}+{C}_{i-2}\right){Q}_{i-1}^{*}\right\}\\ =\frac{1}{{B}_{p+1}+{C}_{p}}\left\{x-\frac{{B}_{p}+{C}_{p-1}}{{C}_{p}}\right\}+\frac{1}{{B}_{p+1}+{C}_{p}}\left\{-{R}_{p}-\underset{i=0}{\overset{p}{\sum }}\left({B}_{i-1}+{C}_{i-2}\right){Q}_{i-1}^{*}\right\}\\ ={Q}_{p+1}+{R}_{p+1},\end{array}$

where

$\begin{array}{l}{Q}_{p+1}:=\frac{1}{{B}_{p+1}+{C}_{p}}\left\{x-\frac{{B}_{p}+{C}_{p-1}}{{C}_{p}}\right\}\\ \text{and}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}{R}_{p+1}:=\frac{1}{{B}_{p+1}+{C}_{p}}\left\{-{R}_{p}-\underset{i=0}{\overset{p}{\sum }}\left({B}_{i-1}+{C}_{i-2}\right){Q}_{i-1}^{*}\right\}.\end{array}$ (7)

${Q}_{p+1}\left(x\right)$ is a polynomial of degree 1 and the residual ${R}_{p+1}\in \text{span}\left\{{Q}_{i}^{*};i=0,1,\cdots ,p-1\right\}$.

Equations (6) and (7) suggest that any canonical function is the sum of a polynomial $\in \text{span}\left\{{x}^{i};i\ge 0\right\}$, and a residual $\in \text{span}\left\{{Q}_{i}^{*};i=0,1,\cdots ,p-1\right\}$. This is formulated in the following theorem:

Theorem 2. For all $k\ge 0$, ${Q}_{k+p}^{*}$ can be written as

${Q}_{k+p}^{*}={Q}_{k+p}+{R}_{k+p}$

where $\left\{{Q}_{k};k\ge 0\right\}$ are called polynomials given by the recursion:

${Q}_{k+p}=\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(k\right)}{Q}_{k+i-2}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\right)$ (8)

with ${Q}_{j}=0$ for all $j=0,1,\cdots ,p-1$. And

${R}_{k+p}=\underset{j=0}{\overset{p-1}{\sum }}\text{ }{r}_{k+p}^{\left(j\right)}{Q}_{j}^{*}$

where $\left\{{r}_{k}^{\left(j\right)};k\ge 0\right\}$ is a sequence of constants given by the recursion

${r}_{k+p}^{\left(j\right)}=-\frac{1}{{f}_{p}^{\left(k\right)}}\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(k\right)}{r}_{k+i-2}^{\left(j\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\text{\hspace{0.17em}}\text{ }\text{and}\text{\hspace{0.17em}}\text{ }j=0,1,\cdots ,p-1\right)$ (9)

with ${r}_{i}^{\left(j\right)}={\delta }_{ij}$ for all i and $j=0,1,\cdots ,p-1$.

Proof. Equations (6)-(7) show that this theorem holds for $k=0$ and $k=1$. The use of an induction argument allows to prove it for all k. To this end let us assume that

${Q}_{k+j}^{*}={Q}_{k+j}+{R}_{k+j},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(j=0,1,2,\cdots ,p-1\right)$ (10)

with

${R}_{j}=\underset{i=0}{\overset{p-1}{\sum }}\text{ }{r}_{j}^{\left(i\right)}{Q}_{i}^{*}$

and let us take $j=p$. Combining (10) with (3) gives:

$\begin{array}{c}{Q}_{k+p}^{*}=\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}\left({Q}_{k+j-2}+{R}_{k+j-2}\right)\right\}\\ =\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}{Q}_{k+j-2}\right\}+\frac{1}{{f}_{p}^{\left(k\right)}}\left\{-\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}{R}_{k+j-2}\right\}\\ =\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}{Q}_{k+j-2}\right\}+\frac{1}{{f}_{p}^{\left(k\right)}}\left\{-\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}\left(\underset{i=0}{\overset{p-1}{\sum }}\text{ }{r}_{k+j-2}^{\left(i\right)}{Q}_{i}^{*}\right)\right\}\\ =\frac{1}{{f}_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}{Q}_{k+j-2}\right\}+\frac{1}{{f}_{p}^{\left(k\right)}}\left\{-\underset{i=0}{\overset{p-1}{\sum }}\left[\underset{j=0}{\overset{p+1}{\sum }}\text{ }{f}_{j-2}^{\left(k\right)}\left({r}_{k+j-2}^{\left(i\right)}\right)\right]{Q}_{i}^{*}\right\}\\ \equiv {Q}_{k+p}+{R}_{k+p}\end{array}$

as required.

Corollary 1. Suppose that the above assumptions and notation hold true.

1) If $p=1$, then $\left\{{Q}_{k};k\ge 0\right\}$ are generated by the recursion:

${Q}_{0}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{Q}_{k+1}=\frac{1}{{f}_{1}^{\left(k\right)}}\left\{{x}^{k}-{f}_{0}^{\left(k\right)}{Q}_{k}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\right).$

And

${R}_{k+1}={r}_{k+1}^{\left(0\right)}{Q}_{0}^{*},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\right)$

where

${r}_{0}^{\left(0\right)}=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{k+1}^{\left(0\right)}=-\frac{{f}_{0}^{\left(k\right)}}{{f}_{1}^{\left(k\right)}}{r}_{k}^{\left(0\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\right)$

2) If $p=2$, then $\left\{{Q}_{k};k\ge 0\right\}$ are generated by the recursion:

${Q}_{0}={Q}_{1}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{Q}_{k+2}=\frac{1}{{f}_{2}^{\left(k\right)}}\left\{{x}^{k}-{f}_{0}^{\left(k\right)}{Q}_{k}-{f}_{1}^{\left(k\right)}{Q}_{k+1}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\right).$

And

${R}_{k+2}={r}_{k+2}^{\left(0\right)}{Q}_{0}^{*}+{r}_{k+2}^{\left(1\right)}{Q}_{1}^{*}$

where $\left\{{r}_{k}^{\left(0\right)},{r}_{k}^{\left(1\right)};k\ge 0\right\}$ are constants given by the recursion

$\begin{array}{l}{r}_{0}^{\left(0\right)}={r}_{1}^{\left(1\right)}=1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{0}^{\left(1\right)}={r}_{1}^{\left(0\right)}=0,\\ {r}_{k+2}^{\left(j\right)}=-\frac{1}{{f}_{1}^{\left(k\right)}}\left({f}_{0}^{\left(k\right)}{r}_{k}^{\left(j\right)}+{f}_{0}^{\left(k\right)}{r}_{k+1}^{\left(j\right)}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(k\ge 0\text{\hspace{0.17em}}\text{ }\text{and}\text{\hspace{0.17em}}\text{ }j=0,1\right).\end{array}$

3. Solutions of Quasi-Exactly Solvable Models in Terms of {Qn}

Having obtained Algorithm (8)-(9) that generates the canonical polynomials that fit to our context, we can now formulate the main result of this paper:

Theorem 3. Let $A\left(x\right)=\underset{i=0}{\overset{p+2}{\sum }}{A}_{i}{x}^{i}$ and $B\left(x\right)=\underset{i=0}{\overset{p+1}{\sum }}{B}_{i}{x}^{i}$ be two given polynomials of degree $p+2$ and $p+1$ respectively, and suppose that $C\left(x\right)=\underset{i=0}{\overset{p}{\sum }}{C}_{i}{x}^{i}$ is an unknown polynomial of degree $p\ge 0$. Let $n\ge 0$. If $\left\{{C}_{0},{C}_{1},\cdots ,{C}_{p}\right\}$ satisfy the following system of algebraic equations

${C}_{p}=-n\left(n-1\right){A}_{p+2}-n{B}_{p+1}$ (11)

$\underset{i=0}{\overset{p+1}{\sum }}\left(n\left(n-1\right){A}_{i}+n{B}_{i-1}+{C}_{i-2}\right){r}_{n+i-2}^{\left(l\right)}=0,\text{ }\left(l=0,1,2,\cdots ,p-1\right)$ (12)

where $\left\{{r}_{i}^{\left(l\right)}\right\}$ are given by (9), then the polynomial

$y\left(x\right)={x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\left(n\left(n-1\right){A}_{i}+n{B}_{i-1}+{C}_{i-2}\right){Q}_{n+i-2}$ (13)

is an exact solution for the differential Equation (2), where $\left\{{Q}_{k}\right\}$ are given by (8).

Proof. Let $n\ge 0$. If (11) holds true, then ${f}_{p}^{\left(n\right)}:=n\left(n-1\right){A}_{p+2}+n{B}_{p+1}+{C}_{p}=0$ and therefore setting $k=n$ in (5) gives

$D{x}^{n}={f}_{p}^{\left(n\right)}{x}^{n+p}+D\left\{\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(n\right)}{Q}_{n+i-2}^{*}\right\}=D\left\{\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i-2}^{\left(n\right)}{Q}_{n+i-2}^{*}\right\}.$

This, in turn, implies that

$D\left\{{x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\left[n\left(n-1\right){A}_{i}+n{B}_{i-1}+{C}_{i-2}\right]{Q}_{n+i-2}^{*}\right\}=0.$

That is,

$y\left(x\right)={x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}{Q}_{n+i-2}^{*}\equiv {x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}\left({Q}_{n+i-2}+{R}_{n+i-2}\right)$ (14)

is an exact solution for the differential Equation (2). For Expression (14) to be a polynomial, it must be independent of the residual terms. To this end, we proceed as follows:

$\begin{array}{c}y\left(x\right)={x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}\left({Q}_{n+i-2}+{R}_{n+i-2}\right)\\ ={x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}\left({Q}_{n+i-2}+\underset{l=0}{\overset{p-1}{\sum }}\text{ }{r}_{n+i-2}^{\left(l\right)}{Q}_{l}^{*}\right)\\ ={x}^{n}-\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}{Q}_{n+i-2}+\underset{l=0}{\overset{p-1}{\sum }}\left(\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}{r}_{n+i-2}^{\left(l\right)}\right){Q}_{l}^{*}\end{array}$

Now, for y(x) to be independent of the p undefined elements $\left\{{Q}_{0}^{*},{Q}_{1}^{*},\cdots ,{Q}_{p-1}^{*}\right\}$, the coefficients of the latter must be set equal to zero leading to the following algebraic system that consists of p equations:

$\underset{i=0}{\overset{p+1}{\sum }}\text{ }{f}_{i}^{\left(n\right)}{r}_{n+i-2}^{\left(l\right)}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(l=0,1,2,\cdots ,p-1\right).$

Explicitly,

$\underset{i=0}{\overset{p+1}{\sum }}\left[n\left(n-1\right){A}_{i}+n{B}_{i-1}+{C}_{i-2}\right]{r}_{n+i-2}^{\left(l\right)}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(l=0,1,2,\cdots ,p-1\right)$

with the p unknowns $\left\{{C}_{0},{C}_{1},{C}_{2},\cdots ,{C}_{p-1}\right\}$. This completes the proof of the theorem. $\square$

For illustration let us consider particular cases:

Case $p=1$. Then we have

$A\left(x\right)=\underset{i=0}{\overset{3}{\sum }}{A}_{i}{x}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}B\left(x\right)=\underset{i=0}{\overset{2}{\sum }}{B}_{i}{x}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}C\left(x\right)=\underset{i=0}{\overset{1}{\sum }}{C}_{i}{x}^{i}$

where $\left\{{C}_{0},{C}_{1}\right\}$ satisfy the algebraic equation

${C}_{1}=n\left(1-n\right){A}_{3}-n{B}_{2}$

$n\left(n-1\right)\left({A}_{0}{r}_{n-2}^{\left(0\right)}+{A}_{1}{r}_{n-1}^{\left(0\right)}+{A}_{2}{r}_{n}^{\left(0\right)}\right)+n\left({B}_{0}{r}_{n-1}^{\left(0\right)}+n{B}_{2}{r}_{n}^{\left(0\right)}\right)+{C}_{0}{r}_{n}^{\left(0\right)}=0$

where $\left\{{r}_{i}^{\left(0\right)}\right\}$ are given by,

$\begin{array}{l}\left(k\left(k-1\right){A}_{3}+k{B}_{2}+{C}_{1}\right){r}_{k+1}^{\left(0\right)}\\ =-k\left(k-1\right){A}_{0}{r}_{k-2}^{\left(0\right)}-\left(k\left(k-1\right){A}_{1}+k{B}_{0}\right){r}_{k-1}^{\left(0\right)}-\left(k\left(k-1\right){A}_{2}+k{B}_{1}+{C}_{0}\right){r}_{k}^{\left(0\right)}\end{array}$

with ${r}_{0}^{\left(0\right)}=1$. And

$y\left(x\right)={x}^{n}-\underset{i=0}{\overset{2}{\sum }}\left(n\left(n-1\right){A}_{i}+n{B}_{i-1}+{C}_{i-2}\right){Q}_{n+i-2}$

where $\left\{{Q}_{k};0\le k\le n\right\}$ are given by the recursion:

$\begin{array}{l}\left(k\left(k-1\right){A}_{3}+k{B}_{2}+{C}_{1}\right){Q}_{k+1}\\ ={x}^{k}-k\left(k-1\right){A}_{0}{Q}_{k-2}-\left(k\left(k-1\right){A}_{1}+k{B}_{0}\right){Q}_{k-1}-\left(k\left(k-1\right){A}_{2}+k{B}_{1}+{C}_{0}\right){Q}_{k}\end{array}$

for $k=0,1,2,\cdots$, with ${Q}_{0}=0$.

Case $p=2$. In this case,

$A\left(x\right)=\underset{i=0}{\overset{4}{\sum }}{A}_{i}{x}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}B\left(x\right)=\underset{i=0}{\overset{3}{\sum }}{B}_{i}{x}^{i},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}C\left(x\right)=\underset{i=0}{\overset{2}{\sum }}{C}_{i}{x}^{i}$

where ${C}_{2}=n\left(1-n\right){A}_{4}-n{B}_{3}$ and $\left\{{C}_{0},{C}_{1}\right\}$ satisfy the $2×2$ algebraic system

$\begin{array}{l}\underset{i=0}{\overset{3}{\sum }}\text{ }n\left(n-1\right){A}_{i}{r}_{n+i-2}^{\left(0\right)}+\underset{i=1}{\overset{3}{\sum }}\text{ }n{B}_{i-1}{r}_{n+i-2}^{\left(0\right)}+\underset{i=2}{\overset{3}{\sum }}\text{ }{C}_{i-2}{r}_{n+i-2}^{\left(0\right)}=0,\\ \underset{i=0}{\overset{3}{\sum }}\text{ }n\left(n-1\right){A}_{i}{r}_{n+i-2}^{\left(1\right)}+\underset{i=1}{\overset{3}{\sum }}\text{ }n{B}_{i-1}{r}_{n+i-2}^{\left(1\right)}+\underset{i=2}{\overset{3}{\sum }}\text{ }{C}_{i-2}{r}_{n+i-2}^{\left(1\right)}=0,\end{array}$

where $\left\{{r}_{k}^{\left(0\right)},{r}_{k}^{\left(1\right)};k\ge 0\right\}$ are given by the recursion

$\begin{array}{l}\left(k\left(k-1\right){A}_{4}+k{B}_{3}+{C}_{2}\right){r}_{k+2}^{\left(j\right)}\\ =-\left({C}_{1}+k{B}_{2}+k\left(k-1\right){A}_{3}\right){r}_{k+1}^{\left(j\right)}-\left({C}_{0}+k{B}_{1}+k\left(k-1\right){A}_{2}\right){r}_{k}^{\left(j\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }-\left(k{B}_{0}+k\left(k-1\right){A}_{1}\right){r}_{k-1}^{\left(j\right)}-k\left(k-1\right){A}_{0}{r}_{k-2}^{\left(j\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}j=0,1\end{array}$

And the exact solution $y\left(x\right)$ is given as:

$\begin{array}{c}y\left(x\right)={x}^{n}-n\left(n-1\right){A}_{0}{Q}_{n-2}-\left[n\left(n-1\right){A}_{1}+n{B}_{0}\right]{Q}_{n-1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left[n\left(n-1\right){A}_{2}+n{B}_{1}+{C}_{0}\right]{Q}_{n}-\left[n\left(n-1\right){A}_{3}+n{B}_{2}+{C}_{1}\right]{Q}_{n+1}\end{array}$

where $\left\{{Q}_{k};0\le k\le n+1\right\}$ are given by the recursion:

$\begin{array}{l}\left(k\left(k-1\right){A}_{4}+k{B}_{3}+{C}_{2}\right){Q}_{k+2}\\ ={x}^{k}-\left[{C}_{1}+k{B}_{2}+k\left(k-1\right){A}_{3}\right]{Q}_{k+1}-\left[{C}_{0}+k{B}_{1}+k\left(k-1\right){A}_{2}\right]{Q}_{k}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }-\left[k{B}_{0}+k\left(k-1\right){A}_{1}\right]{Q}_{k-1}-k\left(k-1\right){a}_{0}{Q}_{k-2}\end{array}$

4. Two Coulombically Repelling Electrons on a Sphere

This is a system governed by the ODE

$\left(\frac{{s}^{2}}{4{R}^{2}}-1\right)\frac{{\text{d}}^{2}U}{\text{d}{s}^{2}}+\left(\frac{\delta s}{4{R}^{2}}-\frac{1}{\gamma s}\right)\frac{\text{d}U}{\text{d}s}+\left(\frac{1}{s}-E\right)U\left(s\right)=0$ (15)

where s denotes the inter-electronic distance between two electrons interacting via a Coulomb potential constrained to remain on the surface of a D-dimensional sphere ( $D\ge 2$ ) of radius R, $\delta =2D-1$, $\gamma =\frac{1}{D-1}$, E is the unknown energy, and $U\left(s\right)$ stands for the unknown inter-electronic wave function (see [8] [9] ).

Setting $x=\frac{s}{2R}$ and $U\left(s\right)=U\left(2Rx\right)\equiv Y\left(x\right)$, then Equation (15) can be written in the form (2) with $p=1$ :

$\left({x}^{3}-x\right)\frac{{\text{d}}^{2}Y}{\text{d}{x}^{2}}+\left(\delta {x}^{2}-\frac{1}{\gamma }\right)\frac{\text{d}Y}{\text{d}x}+\left({C}_{0}+{C}_{1}x\right)Y\left(x\right)=0$ (16)

Let us apply the Algorithm (11)-(12) developed in the previous section to this example. Here

$\begin{array}{l}A\left(x\right)={x}^{3}-x,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{A}_{0}={A}_{2}=0,\text{\hspace{0.17em}}{A}_{1}=-{A}_{3}=-1,\\ B\left(x\right)=\delta {x}^{2}-\frac{1}{\gamma },\text{\hspace{0.17em}}\text{\hspace{0.17em}}{B}_{0}=-\frac{1}{\gamma },\text{\hspace{0.17em}}{B}_{1}=0,\text{\hspace{0.17em}}{B}_{2}=\delta ,\\ C\left(x\right)={C}_{0}+{C}_{1}x,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{C}_{0}=2R,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{C}_{1}=-4{R}^{2}E.\end{array}$ (17)

Let us determine ${C}_{0}$ and ${C}_{1}$, and thereafter obtain R and E. The implementation of Algorithm (11)-(12) yields the following algebraic system that consists of two nonlinear equations with unknowns $\left\{{C}_{0},{C}_{1}\right\}$ :

${C}_{1}=n\left(1-n\right)-n\delta$ (18)

$-n\left(n-1\right){r}_{n-1}^{\left(0\right)}+n\left(\left(-1/\gamma \right){r}_{n-1}^{\left(0\right)}+n\delta {r}_{n}^{\left(0\right)}\right)+{C}_{0}{r}_{n}^{\left(0\right)}=0$ (19)

where $\left\{{r}_{i}^{\left(0\right)}\right\}$ are obtained by,

$\left(k\left(k-1\right)+k\delta +{C}_{1}\right){r}_{k+1}^{\left(0\right)}=\left(k\left(k-1\right)+k/\gamma \right){r}_{k-1}^{\left(0\right)}-{C}_{0}{r}_{k}^{\left(0\right)}$

which is equivalent to the explicit recursion

${r}_{k+1}^{\left(0\right)}=\frac{\left(k\left(k-1\right)+k/\gamma \right){r}_{k-1}^{\left(0\right)}-{C}_{0}{r}_{k}^{\left(0\right)}}{\left(k\left(k-1\right)+k\delta +{C}_{1}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{with}\text{\hspace{0.17em}}{r}_{0}^{\left(0\right)}=1$

Then the desired solution $Y\left(x\right)$ follows from Equation (13):

$Y\left(x\right)={x}^{n}-\underset{i=0}{\overset{2}{\sum }}\left(n\left(n-1\right){A}_{i}+n{B}_{i-1}+{C}_{i-2}\right){Q}_{n+i-2}$

where $\left\{{Q}_{k};0\le k\le n\right\}$ are given by the recursion:

$\left(k\left(k-1\right)+k\delta +{C}_{1}\right){Q}_{k+1}={x}^{k}+\left(k\left(k-1\right)+k/\gamma \right){Q}_{k-1}-{C}_{0}{Q}_{k}$

or, equivalently,

${Q}_{k+1}=\frac{{x}^{k}+\left(k\left(k-1\right)+k/\gamma \right){Q}_{k-1}-{C}_{0}{Q}_{k}}{\left(k\left(k-1\right)+k\delta +{C}_{1}\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{with}\text{\hspace{0.17em}}\text{ }{Q}_{0}=0.$

To illustrate let us consider different values for n.

Case $n=2$ : In this case ${C}_{1}=-4D$ and the values of ${C}_{0}$ are the roots of the equation:

$-{c}_{0}\left({c}_{0}^{2}+2\left(1-4D\right)D\right)=0$

which are $\left\{0,±\sqrt{2D\left(4D-1\right)}\right\}$.

When ${C}_{0}=\sqrt{2D\left(4D-1\right)}$ we obtain:

$E=\frac{2}{4D-1}$

$R=\sqrt{\frac{4{D}^{2}-D}{2}}$

$C\left(x\right)=\sqrt{2D\left(4D-1\right)}-4Dx$

$Y\left(x\right)=\frac{D-1}{2D+1}+\frac{\sqrt{2D\left(4D-1\right)}}{2D+1}x+{x}^{2}$

$\Psi \left(s\right)=Y\left(\frac{s}{2R}\right)=\frac{D-1}{2D+1}+\frac{1}{2D+1}s-\frac{1}{2D-8{D}^{2}}{s}^{2}.$

The value ${C}_{0}=0$ implies that $R=0$ and if ${C}_{0}=-\sqrt{2D\left(4D-1\right)}$ lead to $R<0$. For these reasons, they are discarded.

Case $n=3$ : In this case ${C}_{1}=-3-6D$ and the values of ${C}_{0}$ are the roots of the equations

$-{c}_{0}^{4}+{c}_{0}^{2}\left(20{D}^{2}+20D+6\right)+9\left(-4{D}^{4}-8{D}^{3}+{D}^{2}+8D+3=0\right)$

which are

$\left\{±\sqrt{10{D}^{2}+10D+3±\Delta }\right\}$ where $\Delta =\sqrt{64{D}^{4}+128{D}^{3}+169{D}^{2}+132D+36}$

When ${C}_{0}=\sqrt{10{D}^{2}+10D+3±\Delta }$ we obtain:

$E=\frac{1}{2}\sqrt{10{D}^{2}+10D+3±\Delta }$

$R=\frac{6D+3}{10{D}^{2}+10D+3±\Delta }$

$C\left(x\right)=\left(\sqrt{10{D}^{2}+10D+3±\Delta }\right)-\left(6D+3\right)x$

$\begin{array}{c}Y\left(x\right)=\left[\frac{\left(-4{D}^{2}±\Delta -13D-6\right)\sqrt{10{D}^{2}±\Delta +10D+3}}{12\left(D+1\right)\left(2D+1\right)\left(2D+3\right)}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left[\frac{4{D}^{2}±\Delta -5D-6}{4\left(D+1\right)\left(2D+3\right)}\right]x+\frac{\sqrt{10{D}^{2}±\Delta +10D+3}}{2D+3}{x}^{2}+{x}^{3}\end{array}$

$\Psi \left(s\right)=Y\left(\frac{s}{2R}\right).$

The values ${C}_{0}=-\sqrt{10{D}^{2}+10D+3±\Delta }$ lead to $R<0$ and therefore they are discarded.

Case $n=4$ : In this case ${C}_{1}=-8-8D$ and the values of ${C}_{0}$ are the roots $\xi$ of the equation

$-\xi \left({\xi }^{4}-2{\xi }^{2}\left(20{D}^{2}+45D+28\right)+8\left(32{D}^{4}+144{D}^{3}+185{D}^{2}+18D-64\right)\right)=0$

which are

$\left\{0,\text{\hspace{0.17em}}±\sqrt{20{D}^{2}+45D+28±3\Delta }\right\}$ where

$\Delta =\sqrt{16{D}^{4}+72{D}^{3}+185{D}^{2}+264D+144}$

When ${C}_{0}=\sqrt{20{D}^{2}+45D+28±3\Delta }$ we obtain:

$E=\frac{1}{2}\sqrt{20{D}^{2}+45D+28+3\Delta }$

$R=\frac{8\left(D+1\right)}{20{D}^{2}+3\Delta +45D+28}$

$C\left(x\right)=\left(\sqrt{20{D}^{2}+45D+28±3\Delta }\right)-8\left(D+1\right)x$

$\begin{array}{c}Y\left(x\right)=\left[\frac{\left(D-1\right)\left(\Delta -9D-12\right)}{4\left(D+2\right)\left(2D+3\right)\left(2D+5\right)}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left[\frac{\left(\Delta -9D-12\right)\sqrt{20{D}^{2}+3\Delta +45D+28}}{4\left(D+2\right)\left(2D+3\right)\left(2D+5\right)}\right]x\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left[\frac{3\left(4{D}^{2}+\Delta +3D-4\right)}{4\left(D+2\right)\left(2D+5\right)}\right]{x}^{2}+\left[\frac{\sqrt{20{D}^{2}+3\Delta +45D+28}}{2D+5}\right]{x}^{3}+{x}^{4}\end{array}$

$\Psi \left(s\right)=Y\left(\frac{s}{2R}\right).$

5. Planer Dirac Electron in Coulomb and Magnetic Fields

Such systems are modeled by the covariant Dirac equation (see [10] [11] [12] ),

$i{\gamma }^{\mu }\left({\partial }_{\mu }+ie{A}_{\mu }\right)\Psi \left(t,r\right)={m}_{e}\Psi \left(t,r\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(i=\sqrt{-1}\right)$

which can be transformed to the following second order ODE

$\begin{array}{l}r\left(r+{r}_{0}\right){f}^{″}\left(r\right)+\left(-Be{r}^{3}-Be{r}_{0}{r}^{2}+2\xi r+2L{r}_{0}+{r}_{0}\right){f}^{\prime }\left(r\right)\\ +\left({c}_{2}{r}^{2}+{c}_{1}r+{c}_{0}\right)f\left(r\right)=0\end{array}$ (20)

which is of the form (1) with

$\begin{array}{l}{c}_{2}={E}^{2}-{m}_{e}^{2}-eB\left(\xi +L+1\right){c}_{1}\\ {c}_{1}=2EZ\alpha +\left[{E}^{2}-{m}_{e}^{2}-eB\left(\xi +L+2\right)\right]{r}_{0}{c}_{0}\\ {c}_{0}=2EZ\alpha {r}_{0}+L-\xi \end{array}$

where $L=\mathcal{l}+\frac{1}{2}$ with $\mathcal{l}$ being integer.

This problem was solved in [13] using Bethe Ansatz approach. We will solve it by implementing Algorithm (11)-(12) with $p=2$. The main task is to find expressions for E, Z, B and the exact solution $f\left(r\right)$ as a polynomial of degree n.

Applying Algorithm (11)-(12) we find:

Case $n=1$ : Let $\beta =eB$. Then ${c}_{2}=\beta =eB$, and ${c}_{1}$ and ${c}_{0}$ are the roots of the system:

Coefficient of Q0: $\frac{{c}_{0}{c}_{1}}{\beta }-2{r}_{0}{s}_{0}-{r}_{0}{c}_{0}-{r}_{0}=0$ (21)

Coefficient of Q1: $\frac{{c}_{1}^{2}}{\beta }-{r}_{0}{c}_{1}-2{s}_{0}-{c}_{0}=0$ (22)

First we solve Equation (21) and obtain

${c}_{0}=-\frac{\beta \left(2{r}_{0}\gamma +{r}_{0}\right)}{\beta {r}_{0}-{c}_{1}}.$

Substitute the latter in (22) we arrive to the equivalent cubic equation:

${\beta }^{2}{r}_{0}+\left(2\beta \gamma -{\beta }^{2}{r}_{0}^{2}\right){c}_{1}+2\beta {r}_{0}{c}_{1}^{2}-{c}_{1}^{3}=0.$

Using a computer algebra program such as Mathematica [14] , we find that this cubic equation has one real root only which is:

$\begin{array}{c}{c}_{1}=-\frac{\sqrt[3]{2{\beta }^{3}{r}_{0}{}^{3}-36{\beta }^{2}{r}_{0}\gamma -27{\beta }^{2}{r}_{0}+\Delta }}{3\sqrt[3]{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\sqrt[3]{2}\left(-{\beta }^{2}{r}_{0}{}^{2}-6\beta \gamma \right)}{3\sqrt[3]{2{\beta }^{3}{r}_{0}{}^{3}-36{\beta }^{2}{r}_{0}\gamma -27{\beta }^{2}{r}_{0}+\Delta }}+\frac{2\beta {r}_{0}}{3}\end{array}$

where

$\Delta =\sqrt{4{\left(-{\beta }^{2}{r}_{0}{}^{2}-6\beta \gamma \right)}^{3}+{\left(2{\beta }^{3}{r}_{0}{}^{3}-36{\beta }^{2}{r}_{0}\gamma -27{\beta }^{2}{r}_{0}\right)}^{2}}.$

Further

$f\left(r\right)={r}_{0}-\frac{{c}_{1}}{\beta }+r$

This implies that

${c}_{2}={E}^{2}-{m}_{e}^{2}-eB\left(\xi +L+1\right)=eB$

$\begin{array}{c}{c}_{1}=2EZ\alpha +\left[{E}^{2}-{m}_{e}^{2}-eB\left(\xi +L+2\right)\right]{r}_{0}\\ =-\frac{\sqrt[3]{2{\beta }^{3}{r}_{0}{}^{3}-36{\beta }^{2}{r}_{0}\gamma -27{\beta }^{2}{r}_{0}+\Delta }}{3\sqrt[3]{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\sqrt[3]{2}\left(-{\beta }^{2}{r}_{0}{}^{2}-6\beta \gamma \right)}{3\sqrt[3]{2{\beta }^{3}{r}_{0}{}^{3}-36{\beta }^{2}{r}_{0}\gamma -27{\beta }^{2}{r}_{0}+\Delta }}+\frac{2\beta {r}_{0}}{3}\end{array}$

${c}_{0}=2EZ\alpha {r}_{0}+L-\xi =-\frac{\beta \left(2{r}_{0}\gamma +{r}_{0}\right)}{\beta {r}_{0}-{c}_{1}}.$

Case $n=2$ : For this value of n we find that ${c}_{2}=-2+2\beta$, and ${c}_{1}$ and ${c}_{0}$ are the roots of the system:

$\begin{array}{l}{c}_{0}\left(-2{\beta }^{2}{r}_{0}{}^{2}+\beta \left(4{s}_{0}+2\right)-8{s}_{0}-4\right)+2\left(\beta -1\right){r}_{0}\left(2{s}_{0}+1\right){c}_{1}\\ \text{ }+3\beta {r}_{0}{c}_{0}{c}_{1}+\left(\beta -2\right){c}_{0}^{2}-{c}_{0}{c}_{1}^{2}=4\left(\beta -1\right)\beta {r}_{0}{}^{2}\left(2{s}_{0}+1\right)\end{array}$

$\begin{array}{l}2{c}_{1}\left(-{\beta }^{2}{r}_{0}{}^{2}+4\beta {s}_{0}+\beta -6{s}_{0}-2\right)+3\beta {r}_{0}{c}_{1}^{2}-4\left(\beta -1\right)\beta {r}_{0}{c}_{0}\\ +\left(3\beta -4\right){c}_{0}{c}_{1}-{c}_{1}^{3}=8\left(\beta -1\right){r}_{0}\left(2\beta {s}_{0}+\beta -2\left({s}_{0}+1\right)\right)\end{array}$

which can be solved using Mathematica.

6. Conclusion

In this paper we proposed an efficient method that allows finding closed expressions for the solution of QES models in terms of the canonical polynomials associated with the given differential operator. Unlike the existing method, our method involves solving a system of nonlinear algebraic equations of which the dimensions depend on p and do not increase with n. Several examples were implemented to testify the efficiency of the proposed method.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

 [1] Lanczos, C. (1956) Applied Analysis. Prentice-Hall, Englewood Cliffs, New Jersey. [2] Ortiz, E.L. (1969) The Tau Method. SIAM Journal on Numerical Analysis, 6, 480-492. https://doi.org/10.1137/0706044 [3] Schäfke, R. and Schmidt, D. (1980) The Connection Problem for General Linear Ordinary Differential Equations at Two Regular Singular Points with Applications in the Theory of Special Functions. SIAM Journal on Numerical Analysis, 11, 848-862. https://doi.org/10.1137/0511076 [4] Gómez-Ullate, D., Kamran, N. and Milson, R. (2009) An Extended Class of Orthogonal Polynomials Defined by a Sturm-Liouville Problem. Journal of Mathematical Analysis and Applications, 359, 352-367. https://doi.org/10.1016/j.jmaa.2009.05.052 [5] Bender, C.M. and Dunne, G.V. (1996) Quasi-Exactly Solvable Systems and Orthogonal Polynomials. Journal of Mathematical Physics, 37, 6-11. https://doi.org/10.1063/1.531373 [6] Sasaki, R., Yang, W.L. and Zhang, Y.Z. (2009) Bethe Ansatz Solutions to Quasi-Exactly Solvable Difference Equations. SIGMA, 5, 16. https://doi.org/10.3842/SIGMA.2009.104 [7] Moroz, A. and Miroshnichenko, A.E. (2018) Constraint Polynomial Approach: An Alternative to the Functional Bethe Ansatz Method? arXiv:1807.11871v1 [8] Loos, P.-F. and Gill, P.M.W. (2009) Two Electrons on Hypersphere: A Quasi-Exactly Solvable Model. Physical Review Letters, 103, 123008. https://doi.org/10.1103/PhysRevLett.103.123008 [9] Loos, P.-F. and Gill, P.M.W. (2010) Excited States of Spherium. Molecular Physics, 108, 2527-2532. https://doi.org/10.1080/00268976.2010.508472 [10] Zhang, Y.Z. (2012) Exact Polynomial Solutions of Second Order Differential Equations and Their Applications. Journal of Physics A: Mathematical and Theoretical, 45, 065206. https://doi.org/10.1088/1751-8113/45/6/065206 [11] Agboola, D. and Zhang, Y.-Z. (2012) Unified Derivation of Exact Solutions for a Class of Quasi-Exactly Solvable Models. Journal of Mathematical Physics, 53, 042101. https://doi.org/10.1063/1.3701833 [12] Hatami, N. and Setare, M.R. (2017) Exact Solutions for a Class of Quasi-Exactly Solvable Models: A Unified Treatment. The European Physical Journal Plus, 132, 311. https://doi.org/10.1140/epjp/i2017-11569-6 [13] Chiang, C.-M. and Ho. C.-L. (2002) Planar Dirac Electron in Coulomb and Magnetic Fields: A Bethe Ansatz Approach. Journal of Mathematical Physics, 43, 43-51. https://doi.org/10.1063/1.1418426 [14] Hastings, C., Mischo, K. and Morrison, M. (2016) Hands-On Start to Wolfram Mathematica: And Programming with the Wolfram Language. 2nd Edition, Wolfram Media.