1 Introduction

In a single dimension, a boundary value problem involves an ordinary differential equation (ODE) and conditions that specify the solution and/or its derivatives at a minimum of two points. The quantity of governing conditions matches the order of the ODE. Boundary value problems (BVPs) may be found in various scientific and technological contexts, including but not limited to boundary layer theory in fluid mechanics, heat power transfer theory, and optimization theory. The important book of Agarwal [1] was devoted to studying theoretically the existence and uniqueness theorems for BVPs. In addition, some applications in which these types of equations appear were also presented. Due to the importance of BVPs, they were numerically treated by different techniques. Some algorithms were proposed to treat the linear and non-linear BVPs. For example, numerical algorithms based on Galerkin methods were proposed in [2, 3] to treat some linear BVPs. In contrast, some other algorithms were handled to treat the non-linear BVPs (see, for example, [4, 5]). Among the other schemes that were proposed to treat various BVPs are the differential transform method in [6], operational methods [7], Sinc-collocation method in [8], Green’s function method in [9], iterative methods [10], and spline methods [11]. Orthogonal polynomials (OPs), a topic linked to numerous crucial areas of analysis, have witnessed significant development in recent years. The fields of DEs and integral equations rely heavily on OPs. Additionally, they provide rather general and illustrative examples of specific circumstances in the theory of orthogonal systems. Among the most utilized polynomials in different applications are the Jacobi polynomials (see, for example, [12, 13]). Four kinds of CPs are special polynomials of Jacobi polynomials. It is well-known that these polynomials have their crucial parts in approximation theory. For some contributions that approximate solutions are expressed in terms of different CPs, see, for example, [14, 15].

The use of spectral approaches dramatically influences the extent of numerical analysis; see [16,17,18] for an illustration of their importance in engineering and fluid dynamics. The spectral method’s fundamental principle is to treat approximations as the sum of a few chosen special functions or polynomials. There are three methods to get numerical solutions to differential and integral equations, and their variants are widely used. The three variations differ in the choice of trial and test functions. Using the underlying conditions of the problem as a guide, the Galerkin technique is conditionally applied by choosing the trial functions that meet those conditions; see, for example, [19,20,21]. A critical difference between the tau and Galerkin approaches is that the tau approach is not restricted in choosing basis functions. The collocation method is widely used because it may be used for every differential or integral equation; see for example [22,23,24].

Due to the great importance of the different OPs in various applications, particularly in numerical analysis, their investigations are vital. Recently, Abd-Elhameed in [25] studied a type of OPs, generalized Jacobi polynomials. This paper derived connections between these kinds of OPs and other polynomials. It is worth noting that some studies were performed on some modified CPs to be capable of treating some types of DEs. In [26], the authors introduced certain modified CPs of the second kind and used them to treat the third-order Emden-Fowler type DEs, while the authors in [27] introduced modified CPs of the third-kind (MCPs3K) and employed them to provide solutions to some singular high-order DEs.

Operational matrices of integrals and derivatives are useful in numerical analysis. These matrices are essential to finding numerical solutions to different DEs. (see, for example, [28, 29]). Obtaining an operational matrix of derivatives of certain basis functions relies on expressing the first derivative of this basis in terms of their originals.

Currently, the article introduces three kinds of MCPs3K that will be utilized to solve three types of singular BVPs using the typical collocation method. A unified approach will be proposed for solving the three singular BVPs. Consequently, this paper contributes to:

  • Introducing three new types of modified CPs.

  • Establishing the expressions of the first-order derivative of the modified polynomials.

  • Constructing three new operational matrices of derivatives of the modified polynomials.

  • Designing a unified approach based on applying the typical collocation algorithm to treat three types of singular high-order BVPs.

  • Investigating the convergence of the three expansions of the modified polynomials.

  • Examining the accuracy of the presented algorithm via presenting some illustrative examples.

To the best of our knowledge, the established operational matrices are novel. This gives us motivation to utilize them to solve some singular differential equations.

The paper will proceed as described below. An account of third-kind CPs is provided in Sect. 2, along with introducing three modified versions of third-kind CPs. Three Galekin operational metrics of the three modified CPs are presented in Sect. 3. To deal with the three distinct kinds of singular high-order BVPs, the standard collocation technique is applied in conjunction with the three Galerkin operational matrices proposed in Sect. 4. The convergence and error analysis of the three modified Chebyshev expansions are examined in Sect. 5 via a unified approach. In Sect. 6, some numerical results are displayed. In the final section, we offer some final thoughts.

2 An Account on Third-Kind of CPs and Three Modified Ones

In this section, we will introduce some characteristics of the MCPs3K. In addition, three kinds of modified CPs of the third-kind will be considered. They will be used as the basis functions in our proposed numerical algorithm.

2.1 An Account on Third-Kind CPs

It is well-known CPs of the third-kind are particular polynomials of the Jacobi polynomials that can be expressed as

$$\begin{aligned} {\mathcal {V}}_{m}(\cos \theta )=\displaystyle \frac{\cos \left( \left( m+\frac{1}{2}\right) \theta \right) }{\cos \left( \tfrac{\theta }{2}\right) },\, \theta \in [0,\pi ]. \end{aligned}$$
(1)

Below is a recurrence relation that can be utilized to construct these polynomials:

$$\begin{aligned} {\mathcal {V}}_{m}(t)=2\, t\,{\mathcal {V}}_{m-1}(t)-{\mathcal {V}}_{m-2}(t),\quad {\mathcal {V}}_{0}(t)=1,\quad {\mathcal {V}}_{1}(t)=2\, t-1,\quad m\ge 2. \end{aligned}$$
(2)

The shifted CPs on [ab] can be defined as

$$\begin{aligned} {\mathcal {V}}^{*}_{m}(t)={\mathcal {V}}_{m}\left( \displaystyle \frac{2t-a-b}{b-a}\right) ,\,\,t\in [a,b]. \end{aligned}$$
(3)

\({\mathcal {V}}^{*}_{m}(t),\ m\ge 0\) are orthogonal on [ab] whose orthogonality relation is

$$\begin{aligned} \int \limits _{a}^{b}{w_{0}}(t)\,{\mathcal {V}}^{*}_{i}(t)\,{\mathcal {V}}^{*}_{m}(t)\, dt= {\left\{ \begin{array}{ll} \displaystyle \frac{\pi (b-a)}{2},&{} i=m,\\ 0,&{}i\not =m, \end{array}\right. } \end{aligned}$$
(4)

where \({w_{0}}(t)=\sqrt{\displaystyle \frac{t-a}{b-t}}\).

2.2 Three Kinds of MCPs3K

This section introduces three kinds of OPs that are considered modifications of CPs. The main reason for introducing these polynomials is that these polynomials will meet the underlying conditions of the singular BVPs under investigation.

For our present purposes, and for every positive integer n with \(n\ge 2\), we introduce the three sets of OPs \(\left\{ \phi _{n,j}(t)\right\} _{j\ge 0}\), \(\left\{ \psi _{n,j}(t)\right\} _{j\ge 0}\) and \(\left\{ \chi _{n,j}(t)\right\} _{j\ge 0}\) defined as

$$\begin{aligned} \phi _{n,j}(t)&=(t-a)^{n-1}(b-t)\, {\mathcal {V}}^{*}_{j}(t),\end{aligned}$$
(5)
$$\begin{aligned} \psi _{n,j}(t)&=(t-a)^n(b-t)^n\, {\mathcal {V}}^{*}_{j}(t),\end{aligned}$$
(6)
$$\begin{aligned} \chi _{n,j}(t)&=(t-a)^{n+1}(b-t)^{n}\, {\mathcal {V}}^{*}_{j}(t). \end{aligned}$$
(7)

Remark 1

The three selected basis functions that are defined in (5), (6), and (7) will meet the boundary conditions of the BVPs that will be considered and treated in the subsequent section.

Remark 2

To give a unified approach for our proposed numerical algorithm, we will denote \(\Upsilon _{n,j}(t)\) to be one of the three polynomials \(\phi _{n,j}(t)\), \(\psi _{n,j}(t)\) and \(\chi _{n,j}(t)\).

Remark 3

According to Remark 2 and the orthogonality relation (4), it can be seen \(\Upsilon _{n,j}(t)\) are orthogonal on [ab] regarding \({w_{n}}(t)=\displaystyle \frac{{w_{0}}(t)}{\Omega ^{2}_{n}(t)}\), with \({w_{0}}(t)=\sqrt{\displaystyle \frac{t-a}{b-t}}\), and

$$\begin{aligned} \Omega _{n}(t)={\left\{ \begin{array}{ll} (t-a)^{n-1}(b-t),&{} \Upsilon _{n,i}(t)\equiv \phi _{n,i}(t),\\ (t-a)^{n}(b-t)^{n},&{}\Upsilon _{n,i}(t)\equiv \psi _{n,i}(t),\\ (t-a)^{n+1}(b-t)^{n},&{}\Upsilon _{n,i}(t)\equiv \chi _{n,i}(t). \end{array}\right. } \end{aligned}$$
(8)

Furthermore, we have the following unified orthogonality relation:

$$\begin{aligned} \int _{a}^{b}{w_{n}}(t)\, \Upsilon _{n,i}(t)\, \Upsilon _{n,j}(t)\, dt={\left\{ \begin{array}{ll} \displaystyle \frac{\pi (b-a)}{2},&{} i=j,\\ 0,&{}i\not =j. \end{array}\right. } \end{aligned}$$
(9)

3 Introducing Three New Operational Matrices of Derivatives

Here, we create three new operational matrices of derivatives for the three MCPs3K introduced in (5), (6), and (7). The formulas for these polynomials’ first derivatives provide a path to these operational matrices. Three theorems will be In this concern, we will state and prove three theorems.

Theorem 1

Let \(\phi _{n,i}(t)\) be the modified CPs that are defined in (5). For any integer i with \(i\ge 1\), the first-derivative \(D \phi _{n,i}(t)\) can be expressed as

$$\begin{aligned} D \phi _{n,i}(t)=\displaystyle \sum _{j=0}^{i-1}{\bar{\vartheta }}_{i,j}(n) \phi _{n,j}(t) +{\bar{\delta }}_{n,i}(t), \end{aligned}$$
(10)

where

$$\begin{aligned} {\bar{\vartheta }}_{i,j}(n)= \left( 1+(-1)^{i+j}\right) {\bar{\vartheta }}^{(1)}_{i,j}(n)+ \left( 1-(-1)^{i+j}\right) {\bar{\vartheta }}^{(2)}_{i,j}(n), \end{aligned}$$
(11)

and

$$\begin{aligned} {\bar{\delta }}_{n,i}(t)=(t-a)^{n-2} \left( a-t+(-1)^i (2 i+1) (n-1) (b-t)\right) , \end{aligned}$$
(12)

with

$$\begin{aligned} {\bar{\vartheta }}^{(1)}_{i,j}(n)=&\dfrac{1}{b-a}(j-i)(2n-5),\\ {\bar{\vartheta }}^{(2)}_{i,j}(n)=&\dfrac{1}{b-a}(1+i+j+2\,i\,n-2\,j\,n). \end{aligned}$$

Proof

Let us first suppose \([a,b]=[-1,1]\). In such case, we have

$$\begin{aligned} \phi _{n,i}(t)=(t+1)^{n-1} (1-t){\mathcal {V}}_{i}(t). \end{aligned}$$
(13)

A process of induction is followed on i. In the case of \(i=1\), it is easy to prove that relation (10) holds. Suppose (10) applies to i and \((i-1)\). The next step is to prove that

$$\begin{aligned} D \phi _{n,i+1}(t)=\displaystyle \sum _{j=0}^{i}{\bar{\vartheta }}_{i+1,j}(n) \phi _{n,j}(t) +{\bar{\delta }}_{n,i+1}(t). \end{aligned}$$
(14)

With reference to (2), it is obvious that

$$\begin{aligned} \phi _{n,i+1}(t)&=2t\,\phi _{n,i}(t)-\phi _{n,i-1}(t),\quad i\ge 1 \end{aligned}$$
(15)
$$\begin{aligned} \phi _{n,0}(t)&=(t+1)^{n-1}(1-t),\quad \phi _{n,1}(t)=(t+1)^{n-1}(1-t){\mathcal {V}}_{1}(t). \end{aligned}$$
(16)

Differentiating (15) with respect to t yields

$$\begin{aligned} D\phi _{n,i+1}(t)=2t\,D\phi _{n,i}(t)+2\phi _{n,i}(t)-D\phi _{n,i-1}(t). \end{aligned}$$
(17)

If the induction hypothesis is applied to \(D\phi _{n,i}(t)\) and \(D\phi _{n,i-1}(t)\) in (17), then we get

$$\begin{aligned} D\phi _{n,i+1}(t)=2t\,\sum _{j=0}^{i-1}{\bar{\vartheta }}_{i,j}(n) \phi _{n,j}(t)-\sum _{j=0}^{i-1}{\bar{\vartheta }}_{i-1,j}(n) \phi _{n,j}(t) +2\phi _{n,i}(t)+2t\,{\bar{\delta }}_{n,i}(t)-{\bar{\delta }}_{n,i-1}(t). \end{aligned}$$
(18)

The substitution of (15) written in the form

$$\begin{aligned} 2t\,\phi _{n,j}(t)=\phi _{n,j+1}(t)-\phi _{n,j-1}(t), \end{aligned}$$
(19)

into relation (18), and using \({\bar{\vartheta }}_{i,i}(n)=0\) for all i, will yield-after some algebraic computations- the following formula

$$\begin{aligned} \begin{aligned} D\phi _{n,i+1}(t)&=2t\,\sum _{j=2}^{i-1}\left( {\bar{\vartheta }}_{i,j-1}(n)-{\bar{\vartheta }}_{i-1,j}(n)+{\bar{\vartheta }}_{i,j+1}(n)\right) \phi _{n,j}(t) \\&+{\bar{\vartheta }}_{i,0}(n) \phi _{n,0}(t)+{\bar{\vartheta }}_{i,i-1}(n)\phi _{n,i}(t)+({\bar{\vartheta }}_{i,1}(n)-{\bar{\vartheta }}_{i-1,0}(n))\phi _{n,0}(t)\\&+({\bar{\vartheta }}_{i,2}(n)-{\bar{\vartheta }}_{i-1,1}(n))\phi _{n,1}(t)\\&+2\phi _{n,i}(t)+2t\,{\bar{\delta }}_{n,i}(t)-{\bar{\delta }}_{n,i-1}(t). \end{aligned} \end{aligned}$$
(20)

Formulae (26) and (27) imply the following relations:

$$\begin{aligned}&{\bar{\vartheta }}_{i,j-1}(n)-{\bar{\vartheta }}_{i-1,j}(n)+{\bar{\vartheta }}_{i,j+1}(n)={\bar{\vartheta }}_{i+1,j}(n),\\&{\bar{\vartheta }}_{i,i-1}(n)={\bar{\vartheta }}_{i+1,i}(n),\\&{\bar{\vartheta }}_{i,1}(n)-{\bar{\vartheta }}_{i-1,0}(n)=1+(-1)^{i},\\&{\bar{\vartheta }}_{i,2}(n)-{\bar{\vartheta }}_{i-1,1}(n)=1+(-1)^{i+1},\\&2t\,{\bar{\delta }}_{n,i}(t)-{\bar{\delta }}_{n,i-1}(t)=(t+1)^{n-2} \left( (-1)^{i+1} (n-1) (t-1) (i (4\, t+2)\right. \\&\left. \quad +2 t-1)-t\,(2 t-1)-2 t+1\right) . \end{aligned}$$

It is not difficult to demonstrate the above relations leads to the following formula:

$$\begin{aligned} \begin{aligned}&{\bar{\vartheta }}_{i,0}(n) \phi _{n,0}(t)+{\bar{\vartheta }}_{i,i-1}(n)\phi _{n,i}(t)+({\bar{\vartheta }}_{i,1}(n)-{\bar{\vartheta }}_{i-1,0}(n))\phi _{n,0}(t)\\&+({\bar{\vartheta }}_{i,2}(n)-{\bar{\vartheta }}_{i-1,1}(n))\phi _{n,1}(t)+2\phi _{n,i}(t)+2t\,{\bar{\delta }}_{n,i}(t)-{\bar{\delta }}_{n,i-1}(t)\\&={\bar{\vartheta }}_{i+1,0}(n)\phi _{n,0}(t)+{\bar{\vartheta }}_{i+1,1}(n)\phi _{n,1}(t)+{\bar{\delta }}_{n,i+1}(t), \end{aligned} \end{aligned}$$
(21)

and accordingly, it can be seen that Eq. (18) turns to the form (14).Now, we can finish the proof by substituting \(\dfrac{2t-a-b}{b-a}\) for t in (14). \(\square\)

The following two theorems give the expressions for the first-order derivative of the two modified Cps that are defined in (6) and (7).

Theorem 2

Let \(\psi _{n,i}(t)\) be the modified CPs that are defined in (5). For any integer i with \(i\ge 1\), \(D\psi _{n,i}(t)\) can be expressed as

$$\begin{aligned} D \psi _{n,i}(t)=\displaystyle \sum _{j=0}^{i-1}{\tilde{\vartheta }}_{i,j}(n) \psi _{n,j}(t) +{\tilde{\delta }}_{n,i}(t), \end{aligned}$$
(22)

where

$$\begin{aligned} {\tilde{\vartheta }}_{i,j}(n)= \left( 1+(-1)^{i+j}\right) {\tilde{\vartheta }}^{(1)}_{i,j}(n)+ \left( 1-(-1)^{i+j}\right) {\tilde{\vartheta }}^{(2)}_{i,j}(n), \end{aligned}$$
(23)

and

$$\begin{aligned} {\tilde{\delta }}_{n,i}(t)=n((t-a) (b-t))^{n-1} \left( a-t+(-1)^i (2 i+1) (b-t)\right) , \end{aligned}$$
(24)

with

$$\begin{aligned} {\tilde{\vartheta }}^{(1)}_{i,j}(n)=&\dfrac{1}{b-a}(i-j),\\ {\tilde{\vartheta }}^{(2)}_{i,j}(n)=&\dfrac{1}{b-a}((4 n+1) (i-j)+(2 j+1)). \end{aligned}$$

Proof

Similarly to the proof of Theorem 1. \(\square\)

Theorem 3

Let \(D \chi _{n,i}(t)\) be the modified CPs that are defined in (5). For any integer i with \(i\ge 1\), \(D \chi _{n,i}(t)\) can be expressed as

$$\begin{aligned} D \chi _{n,i}(t)=\displaystyle \sum _{j=0}^{i-1}{\hat{\vartheta }}_{i,j}(n) \chi _{n,j}(t) +{\hat{\delta }}_{n,i}(t), \end{aligned}$$
(25)

where

$$\begin{aligned} {\hat{\vartheta }}_{i,j}(n)= \left( 1+(-1)^{i+j}\right) {\hat{\vartheta }}^{(1)}_{i,j}(n)+ \left( 1-(-1)^{i+j}\right) {\hat{\vartheta }}^{(2)}_{i,j}(n), \end{aligned}$$
(26)

and

$$\begin{aligned} {\hat{\delta }}_{n,i}(t)= (t-a)^n\,(b-t)^{n-1} \left( n(a-t)+(-1)^i (n+1)(2 i+1) (b-t)\right) , \end{aligned}$$
(27)

with

$$\begin{aligned} {\hat{\vartheta }}^{(1)}_{i,j}(n)=&\dfrac{1}{b-a}(j-i),\\ {\hat{\vartheta }}^{(2)}_{i,j}(n)=&\dfrac{1}{b-a}((4 n+3) (i-j)+2 j+1). \end{aligned}$$

Proof

Similar to the proof of Theorem 1. \(\square\)

Remark 4

Aiming to utilize a unified numerical algorithm for three different BVPs using the three modified CPs in (5), (6), and (7), we will unify the results of Theorems 1,2 and 3 in the following theorem.

Theorem 4

Let \(\Upsilon _{n,i}(t)\) be one of the three modified polynomials that are defined in (5), (6), and (7). \(D \Upsilon _{n,i}(t),\,i\ge 1,\) can be expressed as

$$\begin{aligned} D \Upsilon _{n,i}(t)=\displaystyle \sum _{j=0}^{i-1}\vartheta _{i,j}(n) \Upsilon _{n,j}(t) +\delta _{n,i}(t), \end{aligned}$$
(28)

where

$$\begin{aligned} \vartheta _{i,j}(n)=&{\left\{ \begin{array}{ll} {\bar{\vartheta }}_{{i,j}}(n),&{} \Upsilon _{n,i}(t)\equiv \phi _{n,i}(t), \\ {\tilde{\vartheta }}_{{i,j}}(n),&{} \Upsilon _{n,i}(t)\equiv \psi _{n,i}(t),\\ {\hat{\vartheta }}_{{i,j}}(n),&{} \Upsilon _{n,i}(t)\equiv \chi _{n,i}(t), \end{array}\right. } \end{aligned}$$
(29)

and

$$\begin{aligned} \delta _{n,i}(t)={\left\{ \begin{array}{ll} {\bar{\delta }}_{n,i}(t),&{} \Upsilon _{n,i}(t)\equiv \phi _{n,i}(t),\\ {\tilde{\delta }}_{n,i}(t),&{} \Upsilon _{n,i}(t)\equiv \psi _{n,i}(t),\\ {\hat{\delta }}_{n,i}(t),&{} \Upsilon _{n,i}(t)\equiv \chi _{n,i}(t). \end{array}\right. } \end{aligned}$$
(30)

Proof

Theorem 4 is a direct consequence of Theorems  12 and 3. \(\square\)

3.1 Function Approximation

Now, if consider a function \(u(t)\in \, L^{2}_{w_{n}}[a,b]\) that can be expanded as:

$$\begin{aligned} u(t)=\displaystyle \sum _{i=0}^{\infty }a_i\, \Upsilon _{n,i}(t), \end{aligned}$$
(31)

where

$$\begin{aligned} a_i=\frac{2}{\pi \,(b-a)}\,\int \limits _{a}^{b}w_{n}(t)\,u(t)\, \Upsilon _{n,i}(t)\,dt. \end{aligned}$$
(32)

Let us assume that we can approximate the series in Eq. (31) as

$$\begin{aligned} u(t)\simeq u_{N}(t)=\displaystyle \sum _{i=0}^{N}a_i\, \Upsilon _{n,i}(t)={\varvec{A}}^T\,{\varvec{\Upsilon _{n,N}}}(t), \end{aligned}$$
(33)

where

$$\begin{aligned} {\varvec{A}}^T=\left[ a_0, a_1,\dots ,a_N\right] ,\, {\varvec{\Upsilon _{n,N}}}(t)=[\Upsilon _{n,0}(t), \Upsilon _{n,1}(t),\dots , \Upsilon _{n,N}(t)]^{T}. \end{aligned}$$
(34)

3.2 Three Operational Matrices of Modified CPs

Based on the three derivative formulas of the basis functions considered in (5),(6), and (7) that are unified in Theorem 4, the following corollary presents the operational matrix of the derivative of \({\varvec{\Upsilon _{n,N}}}(t)\).

Corollary 1

In terms of \({\varvec{\Upsilon _{n,N}}}(t),\ \displaystyle \frac{d{\varvec{\Upsilon _{n,N}}}(t)}{dt}\), the following formula holds:

$$\begin{aligned} \displaystyle \frac{d{\varvec{\Upsilon _{n,N}}}(t)}{dt}= {G_{n}}\,{\varvec{\Upsilon _{n,N}}}(t)+{\varvec{\delta _{n}}}(t), \end{aligned}$$
(35)

where \({\varvec{\delta _{n}}}(t)=\left( \delta _{n,0}(t),\delta _{n,1}(t),\dots , \delta _{n,N}(t)\right) ^T and \ G_{n}=\big (g^{(n)}_{ij}\big )_{0\le i,j\le N}\). The entries of \(G_{n}\) have the following form:

$$\begin{aligned}g^{(n)}_{ij}={\left\{ \begin{array}{ll} \vartheta _{i,j}(n)\, &{}\quad i>j,\\ 0,&{}\quad otherwise. \end{array}\right. } \end{aligned}$$

Now, we give the operational matrices of derivatives for the three considered modified CPs as examples. For \(N=6\), the corresponding operational matrices for every \(n\ge 2,\) have the following forms:

  1. (i)

    If \(\Upsilon _{n,i}(t)=\phi _{n,i}(t)\)

    $$\begin{aligned} \hspace{-40pt}G_{n}=\displaystyle \frac{4}{b-a} \begin{pmatrix} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ n+1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 5-2 n &{} n+2 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 3 n+2 &{} 5-2 n &{} n+3 &{} 0 &{} 0 &{} 0 &{} 0 \\ 10-4 n &{} 3 (n+1) &{} 5-2 n &{} n+4 &{} 0 &{} 0 &{} 0 \\ 5 n+3 &{} 10-4 n &{} 3 n+4 &{} 5-2 n &{} n+5 &{} 0 &{} 0 \\ 15-6 n &{} 5 n+4 &{} 10-4 n &{} 3 n+5 &{} 5-2 n &{} n+6 &{} 0 \\ \end{pmatrix} _{7\times 7}, \end{aligned}$$
    (36)
  2. (ii)

    If \(\Upsilon _{n,i}(t)=\psi _{n,i}(t)\)

    $$\begin{aligned} G_{n}=\displaystyle \frac{4}{b-a} \begin{pmatrix} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 2 n+1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 1 &{} 2 (n+1) &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 6 n+2 &{} 1 &{} 2 n+3 &{} 0 &{} 0 &{} 0 &{} 0 \\ 2 &{} 6 n+3 &{} 1 &{} 2 (n+2) &{} 0 &{} 0 &{} 0 \\ 10 n+3 &{} 2 &{} 6 n+4 &{} 1 &{} 2 n+5 &{} 0 &{} 0 \\ 3 &{} 10 n+4 &{} 2 &{} 6 n+5 &{} 1 &{} 2 (n+3) &{} 0 \\ \end{pmatrix} _{7\times 7}, \end{aligned}$$
    (37)
  3. (iii)

    If \(\Upsilon _{n,i}(t)=\chi _{n,i}(t)\)

    $$\begin{aligned} G_{n}=\displaystyle \frac{4}{b-a} \begin{pmatrix} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 2 (n+1) &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ -1 &{} 2 n+3 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 6 n+5 &{} -1 &{} 2 (n+2) &{} 0 &{} 0 &{} 0 &{} 0 \\ -2 &{} 6 (n+1) &{} -1 &{} 2 n+5 &{} 0 &{} 0 &{} 0 \\ 10 n+8 &{} -2 &{} 6 n+7 &{} -1 &{} 2 (n+3) &{} 0 &{} 0 \\ -3 &{} 10 n+9 &{} -2 &{} 6 n+8 &{} -1 &{} 2 n+7 &{} 0 \\ \end{pmatrix}_{7\times 7}. \end{aligned}$$
    (38)

Due to Corollary 1, \(\displaystyle \frac{d^{r}\varvec{\Upsilon _{n,N}}(t)}{dt^{r}}\) can be easily evaluated. This result can be seen in the corollary below.

Corollary 2

For every positive integer r, \(\displaystyle \frac{d^{r}\varvec{\Upsilon _{n,N}}(t)}{dt^{r}}\) can be represented as

$$\begin{aligned} \frac{d^{r}\varvec{\Upsilon _{n,N}}(t)}{dt^{r}}= G^{r}_{n}\varvec{\Upsilon _{n,N}}(t)+\varvec{\eta _{n}}^{(r)}(t), \end{aligned}$$
(39)

with \(\varvec{\eta _{n}}^{(r)}(t)=\sum \limits _{k=0}^{r-1}G^{k}_{n}\,\varvec{\delta _{n}}^{(r-k-1)}(t)\).

Proof

Repeated application of (35) yields (39). \(\square\)

Remark 5

In the upcoming section, we will treat three types of singular BVPs with domain [0, 1]. So, we will employ the three operational matrices established in Sect. 3 for \(a=0,b=1\).

4 Galerkin Operational Matrix Method for Treating Three Types of Singular High-order BVPs

In this section, we show how to use the Galerkin operational matrices of derivatives given in (39) to numerically solve the high-order singular BVPs in the form [30,31,32]

$$\begin{aligned} u^{(n)}(t)+\dfrac{m}{t}\,u^{(n-1)}(t)+Lu(t)=g(t),\quad 0\le t\le 1,\,\,n\ge 2, \end{aligned}$$
(40)

subject to one of the following three types of the BCs:

$$\begin{aligned} u^{(i)}(0)=&\alpha _{i},\,\,i=0,1,\dots ,n-2,\,\,\, \text {and}\quad u(1)=\beta _{0}, \end{aligned}$$
(41)
$$\begin{aligned} u^{(i)}(0)=&\alpha _{i},\,i=0,1,\dots ,\dfrac{n-2}{2},\,\, \text {and}\,\, u^{(i)}(1)=\beta _{i},\,i=0,1,\dots ,\dfrac{n-2}{2}\,\,\, \text {for}\, n\, \text {even}, \end{aligned}$$
(42)
$$\begin{aligned} u^{(i)}(0)=&\alpha _{i} ,\,i=0,1,\dots ,\dfrac{n-1}{2},\,\, \text {and}\,\, u^{(i)}(1)=\beta _{i} ,\,i=0,1,\dots ,\dfrac{n-3}{2}\,\,\, \text {for}\, n \,odd, \end{aligned}$$
(43)

where L is a nonlinear differential operator of an order less than n.

The main idea in the proposed algorithm is to convert (40) subject to the non-homogeneous BCs conditions (41), (42) or (43) into a modified one subject to the homogeneous BCs. For this purpose, we set

$$\begin{aligned} {\tilde{u}}(t)=u(t)-P_{n-1}(t). \end{aligned}$$
(44)

The determination of the polynomials \(P_{n-1}(t)\) the type of the boundary condition governed by the BVP (40). For the BCs (41), \(P_{n-1}(t)\) takes the form

$$\begin{aligned} P_{n-1}(t)=\sum \limits _{k=0}^{n-2}\dfrac{\alpha _{k}}{k!}\,t^{k}+\left( \beta _{0}-\sum \limits _{k=0}^{n-2}\dfrac{\alpha _{k}}{k!}\right) \, t^{n-1}. \end{aligned}$$
(45)

For the BCs (42), and (43), to determine \(P_{n-1}(t)\), it can be assumed that

$$\begin{aligned} P_{n-1}(t)= {\left\{ \begin{array}{ll} \sum \limits _{k=0}^{\frac{n}{2}-1}\dfrac{\alpha _{k}}{k!}\,t^{k} +\sum \limits _{k=\frac{n}{2}}^{n-1}A_{k}\,t^{k}, &{} \text {for the BCs} (42),\\ \sum \limits _{k=0}^{\frac{n-1}{2}}\dfrac{\alpha _{k}}{k!}\,t^{k} +\sum \limits _{k=\frac{n+1}{2}}^{n-1}B_{k}\,t^{k}, &{} \text {for the BCs} (43), \end{array}\right. } \end{aligned}$$
(46)

and the coefficients \(A_{i}\, (i=\frac{n}{2},\dots ,n-1)\) and \(B_{i}\, (i=\frac{n+1}{2},\dots ,n-1)\) can be computed by solving the following two systems of equations:

$$\begin{aligned} \sum \limits _{k=\frac{n}{2}}^{n-1}\dfrac{k!}{(k-i)!}A_{k}= \beta _{i}-\sum \limits _{k=i}^{\frac{n}{2}-1}\dfrac{\alpha _{i}}{(k-i)!},\,\,i=0,\dots ,\frac{n-2}{2}, \end{aligned}$$
(47)

and

$$\begin{aligned} \sum \limits _{k=\frac{n+1}{2}}^{n-1}\dfrac{k!}{(k-i)!}B_{k}= \beta _{i}-\sum \limits _{k=i}^{\frac{n-1}{2}}\dfrac{\alpha _{i}}{(k-i)!},\,\,i=0,\dots ,\frac{n-3}{2}. \end{aligned}$$
(48)

Now, making use of the transformation (44), it can be seen Eq. (40) governed by the on-homogeneous BCs (41), (42) or (43) turns into a modified one governed by the homogeneous BCs. More precisely, the modified equation is

$$\begin{aligned} {\tilde{u}}^{(n)}(t)+\dfrac{m}{t}\,{\tilde{u}}^{(n-1)}(t)+L({\tilde{u}}(t)+P_{n-1}(t))= g(t)-\,\dfrac{C_{n}}{t},\quad 0\le t\le 1, \end{aligned}$$
(49)

subject to one of the following the homogeneous BCs:

$$\begin{aligned} {\tilde{u}}^{(i)}(0)=&0,\,\,i=0,1,\dots ,n-2,\,\, \text {and}\,\,\, {\tilde{u}}(1)=0, \end{aligned}$$
(50)
$$\begin{aligned} {\tilde{u}}^{(i)}(0)=&{\tilde{u}}^{(i)}(1)=0,\,i=0,1,\dots ,\dfrac{n-2}{2},\,\,\, \text {for } n \text { even },\end{aligned}$$
(51)
$$\begin{aligned} {\tilde{u}}^{(i)}(0)=&{\tilde{u}}^{(i)}(1)=0 ,\,i=0,1,\dots ,\dfrac{n-3}{2},\,\, \text {and}\,\, {\tilde{u}}^{(\frac{n-1}{2})}(0)=0 \,\,\, \text { for } n \text { odd }, \end{aligned}$$
(52)

with

$$\begin{aligned} C_{n}= {\left\{ \begin{array}{ll} m\, (n-1)!\, \left( \beta _{0}-\sum \limits _{k=0}^{n-2}\dfrac{\alpha _{k}}{k!}\right) ,&{}\text {for the BCs}\,\, (41),\\ m\,A_{n-1}(n-1)!,&{}\text {for the BCs}\,\, (42),\\ m\,B_{n-1}(n-1)!,&{}\text {for the BCs}\,\, (43). \end{array}\right. } \end{aligned}$$

Now, we will propose a unified collocation approach namely shifted Chebyshev third kind Galerkin collocation operational matrix methods, (SC3GCOMMB1),  (SC3GCOMMB2) and (SC3GCOMMB3) using the three considered basis functions to solve the singular BVPs (49)

$$\begin{aligned} {\left\{ \begin{array}{ll} \phi _{n,j}(t)&{}=t^{n-1}(1-t)\,{\mathcal {V}}^{*}_{j}(t),\\ \psi _{{\frac{n}{2}},j}(t)\,\,\,\,&{}=t^{\frac{n}{2}}(1-t)^{\frac{n}{2}}\, {\mathcal {V}}^{*}_{j}(t),\quad \,\,\,\text {for}\, n\,\text {even}, \\ \chi _{{\frac{n-1}{2}},j}(t)&{}=t^{\frac{n+1}{2}}(1-t)^{{\frac{n-1}{2}}}\, {\mathcal {V}}^{*}_{j}(t), \text {for}\, n\,\text {odd}, \end{array}\right. } \end{aligned}$$
(53)

for the purpose of finding approximate solutions to Eq. (49) directed with the three homogeneous BCs (50), (51) or (52), respectively. Through the use of Eqs. (35) and (39), we are able to determine the high-order derivatives of the approximate solution \(\tilde{{u}}_{N}(t)\). Hence, the residual of Eq. (49) can be written as

$$\begin{aligned} \begin{aligned} R_{n,N}(t) = &\varvec{A}^T\, G^{n}_{q}\, \varvec{\Upsilon _{n,N}}(t)+\varvec{\eta }^{(n)}_{q}(t)+\dfrac{m}{t} (\varvec{A}^T\, G^{n-1}_{q}\, \varvec{\Upsilon _{n,N}}(t)+\varvec{\eta }^{(n-1)}_{q}(t))\\&+L(\varvec{A}^T\,\varvec{\Upsilon _{n,N}}(t)+P_{n-1}(t))-\left( g(t)-\,\dfrac{C_{n}}{t}\right) , \end{aligned} \end{aligned}$$
(54)

where \(q\in \left\{ n,\,\frac{n}{2},\,\frac{n-1}{2}\right\}\), according to the choice of basis (53). The application of shifted Chebyshev third kind collocation operational matrix method to solve numerically (49) subject to one of BCs (50), (51) or (52). The collocation points are selected as the \((N+1)\) zeros of \({\mathcal {V}}^{*}_{N+1}(t)\). Consider them \(t_{i},\, i=0,1,\dots ,N\). Thus

$$\begin{aligned} R_{n,N}(t_{i})=0,\quad i=0,1,...,N. \end{aligned}$$
(55)

Any appropriate algorithm may be used to find the solution \(\tilde{{u}}_{N}(t)\) to the \((N+1)\) nonlinear algebraic equations generated by the system (55).

5 Convergence and Error Analysis of the Three Modified Chebyshev Expansions

Here, we offer the results of our research into the error analysis of the proposed algorithm. Suppose that

$$\begin{aligned} {\tilde{u}}(t)=\Omega _{q}(t)g(t), \end{aligned}$$

where \(q\in \left\{ n,\tfrac{n}{2},\,\tfrac{n-1}{2}\right\}\) according to the considered BCs (50)–(52), respectively, with \(g(t)\in C^{n}[0,1]\le M,\) for all \(x\in [0,1]\), and M is a positive integer.

Theorem 5

Assume that \({\tilde{u}}(t)\in L^n_{w_{q}}[0,1]\) be the exact solution of the differential Eq. (49) subject to one of the BCs (50), (51) or (52). let \({\tilde{u}}(t)\) be expanded as

$$\begin{aligned} {\tilde{u}}(t)=\displaystyle \sum _{i=0}^{\infty }c_{i}\, \Upsilon _{n,i}(t). \end{aligned}$$

The series converges uniformly to \({\tilde{u}}(t)\), and the coefficients \(c_{i}\) satisfy the the following inequality

$$\begin{aligned} |c_{i}|<\frac{M}{2^{n-1}\, (i-n+1)^n},\qquad \forall \ i>n-1. \end{aligned}$$
(56)

Proof

The application of the orthogonality relation (4) yields

$$\begin{aligned} c_{i}=\frac{2}{\pi }\displaystyle \int _{0}^{1} {w_{q}}(t) {\tilde{u}}(t)\Upsilon _{n,i}(t)\,dt =\frac{2}{\pi }\displaystyle \int _{0}^{1} w_{0}(t)g(t){\mathcal {V}}^{*}_{i}(t)\,dt. \end{aligned}$$
(57)

The substitution by: \(2\, t-1=\cos \theta\) in (57) yields

$$\begin{aligned} c_{i}=\frac{2}{\pi }\displaystyle \int _{0}^{\pi } g\left( \frac{1+\cos \theta }{2}\right) \,\cos (\theta /2)\,\cos ((i+1/2))\theta \,d\theta . \end{aligned}$$

The integration by parts n times leads to

$$\begin{aligned} c_{i}=\frac{1}{2^{n-1}\pi }\displaystyle \int _{0}^{\pi } g^{(n)}\left( \frac{1+\cos \theta }{2}\right) \, R_i(\theta )\,d\theta , \end{aligned}$$

where \(R_{i}(\theta )\) is a trigonometric polynomial in \(\cos \theta\)\(\sin \theta\). So, using \(|\cos \theta |\le 1, |\sin \theta |\le 1\) and after some manipulations we get the desired result. \(\square\)

Theorem 6

If we consider an approximate solution of (49) in the form \({\tilde{u}}_{N}(t)=\displaystyle \sum _{i=0}^{N}c_{i}\, \Upsilon _{n,i}(t)\). Under the same conditions stated in Theorem 5, then the global error satisfies the following inequality:

$$\begin{aligned} \Vert y-y_{N}\Vert _{w_{q}}< \displaystyle \frac{M}{2^{n-\frac{3}{2}}\, N^{n-\frac{1}{2}}}. \end{aligned}$$
(58)

Proof

As a consequence of the orthogonality property of \(\Upsilon _{n,i}(t),\, i\ge 0\), we have

$$\begin{aligned} \Vert y-y_{N}\Vert ^2_{w_{q}}=\frac{\pi }{2} \displaystyle \sum _{i=N+1}^{\infty }c^2_{i}, \end{aligned}$$

and accordingly, Theorem 5 leads to

$$\begin{aligned} \Vert y-y_{N}\Vert ^2_{w_{q}}<\frac{M^2\,\pi }{2^{2n-1}} \displaystyle \sum _{i=N+1}^{\infty }\frac{1}{(i-n+1)^{2n}}, \end{aligned}$$

and therefore, the following inequality also holds (see, [33])

$$\begin{aligned} \Vert y-y_{N}\Vert ^2_{w_{q}}<\frac{M^2\,\pi }{2^{2n-1}}\displaystyle \int _N^{\infty }\frac{1}{(t-n+1)^{2n}}\,dt. \end{aligned}$$
(59)

In conclusion, the last inequality suggests that

$$\begin{aligned} \Vert y-y_{N}\Vert _{w_{q}}< \displaystyle \frac{M}{2^{n-\frac{3}{2}}\, N^{n-\frac{1}{2}}}. \end{aligned}$$

Theorem 6 is now proved. \(\square\)

6 Numerical Results

This section aims to demonstrate the usefulness of our suggested collocation method by providing some instances. Additionally, we provide some comparisons with other approaches. Furthermore, in every instance within this section, the mistakes are assessed using the maximum norm, namely \(E_{N}=\max \limits _{0 \le t \le 1}|u(t)-u_{N}(t)|\).

Remark 6

We point out that the numerical examples presented in this section to solve Eq. (40) subject to the given boundary conditions will be obtained using Mathematica 13.3 on a computer system equipped with an Intel(R) Core(TM) i9-10850 CPU operating at 3.60 GHz, featuring 10 cores and 20 logical processors.

Example 1

Consider the following two singular BVPs:

$$\begin{aligned} \begin{aligned}&u^{(5)}(t)+\dfrac{1}{t}u^{(4)}(t)+2\,u^{(2)}(t)+3\,u^{(1)}(t)-(u(t))^{2}=-1+1080\\&\quad t+60\,t^4+18\,t^5-2\,t^6-t^{12}, \end{aligned} \end{aligned}$$
(60)

subject to the BCs:

$$\begin{aligned} u(0)=1,\,\,u^{(i)}(0)=0,\,i=1,2,3,\,\,u(1)=2, \end{aligned}$$
(61)

or

$$\begin{aligned} u(0)=1,\,u^{(1)}(0)=0,\,u^{(2)}(0)=0,\,u(1)=2,\,u^{(1)}(1)=6. \end{aligned}$$
(62)

The exact solution of (60) is: \(u(t)=1+t^{6}.\) Our unified collocation algorithm is applied, and the approximate solution is given by

$$\begin{aligned} {\tilde{u}}_{N}(t)=\displaystyle \sum _{i=0}^{N}c_{i}\, \Upsilon _{n,i}(t), N=1,2,\dots , \end{aligned}$$

where

$$\begin{aligned} \Upsilon _{i}(t)&={\left\{ \begin{array}{ll} t^{4}(1-t)\, {\mathcal {V}}^{*}_{j}(t),&{} \text {for the BCs}\, (61), \\ t^{3}(1-t)^2\, {\mathcal {V}}^{*}_{j}(t),&{}\text {for the BCs}\, (62), \end{array}\right. }\end{aligned}$$
(63)
$$\begin{aligned} c_{0}&={\left\{ \begin{array}{ll} -\frac{7}{4},&{} \Upsilon _{n,i}\equiv t^{4}(1-t)\, {\mathcal {V}}^{*}_{i}(t),\\ \frac{11}{4},&{} \Upsilon _{n,i}\equiv t^{3}(1-t)^2\, {\mathcal {V}}^{*}_{i}(t), \end{array}\right. }\quad c_{1}={\left\{ \begin{array}{ll} -\frac{1}{4},&{} \Upsilon _{n,i}\equiv t^{4}(1-t)\, {\mathcal {V}}^{*}_{i}(t),\\ \frac{1}{4},&{} \Upsilon _{n,i}\equiv t^{3}(1-t)^2\, {\mathcal {V}}^{*}_{i}(t), \end{array}\right. } \end{aligned}$$
(64)

and \(c_{i}=0, i=2,3,\dots ,N\). In view of (44), we obtain \(u_{N}(t)=1+t^6\), which is the exact solution.

Example 2

Consider the singular BVP: [30, 31, 34, 35]

$$\begin{aligned} u^{(3)}(t)-\dfrac{2}{t}u^{(2)}(t)-u(t)-u^{2}(t)=-6\,e^{t}+6\,t\,e^{t}+7t^{2}e^{t}-t^{6}\,e^{2t},\,\,\,\,\,0\le t\le 1, \end{aligned}$$
(65)

subject to the BCs:

$$\begin{aligned} u(0)=0,\,\,u^{(1)}(0)=0,\,\,u(1)=e. \end{aligned}$$

The exact solution is \(u(t)=t^{3}e^{t}\).

We apply our suggested method SC3GCOMMB1 for various values of N. Table 1 shows the maximum errors that occur accompanied by the CPU time (in seconds), and Table 2 compares the results of various built techniques. Furthermore, Fig. 1 presents logerror for various N to show the stability of solutions of Example 2. The application of our proposed method SC3GCOMMB1.

Table 1 Maximum absolute error and CPU time (seconds) of Example 2
Table 2 Comparison between different methods of Example 2
Fig. 1
figure 1

Absolute computational error and LogError of Example 2

Example 3

Consider the following singular BVPs:

$$\begin{aligned}&u^{(4)}(t)+\dfrac{4}{t}u^{(3)}(t)-15u^{5}(t)(3-7t^{2}u^{2}(t))(1-t^{2}u^{2}(t))=g(t),\,\,\,\,0\le t\le 1, \end{aligned}$$
(66)
$$\begin{aligned}&u^{(4)}(t)+\dfrac{4}{t}u^{(3)}(t)-15u^{5}(t)(3-7t^{2}u^{2}(t))(1-t^{2}u^{2}(t))=0,\,\,\,\,0\le t\le 1, \end{aligned}$$
(67)

subject to the BCs:

$$\begin{aligned} u(0)=\dfrac{1}{2},\,u^{(1)}(0)=0,\,u(1)=\dfrac{1}{\sqrt{5}},\,u^{(1)}(1)=-\dfrac{1}{5\sqrt{5}}, \end{aligned}$$

where g(t) is chosen such that the exact solution of (66) is

$$\begin{aligned}u(t)=\dfrac{1}{2}+\left( \dfrac{16}{5\sqrt{5}}-\dfrac{1}{2}\right) t^{2}-\left( 3+\dfrac{11}{5\sqrt{5}}\right) t^{3}+6t^{4}-4t^{5}+t^{6},\end{aligned}$$

and the exact solution of (67) is \(u(t)=\dfrac{1}{\sqrt{4+t^2}}\).

The application of proposed method SC3GCOMMB2 to solve (66) gives

$$\begin{aligned} {\tilde{u}}_{N}(t)=\displaystyle \sum _{i=0}^{N}c_{i}\, \Upsilon _{n,i}(t),\, N=2,3,\dots , \end{aligned}$$

where \(\Upsilon _{i}(t)= t^{2}(1-t)^2\, V^{*}_{j}(t)\). In such case, we have

$$\begin{aligned} c_{0}=\dfrac{125}{1000},c_{1}=\dfrac{-1875}{10000},c_{2}=\dfrac{625}{10000}\,\, and\,\, c_{i}=0,i=3,4,\dots ,N. \end{aligned}$$

The greatest errors obtained by solving (67) using our suggested method SC3GCOMMB2 with various values of N are shown in Table 3, accompanied by the CPU (in seconds), and a comparison of the various techniques is shown in Table  4. Furthermore, Fig. 2a shows the computational error norm with \(N = 13\), and Fig. 2b shows the approximate and precise solutions.

Table 3 Maximum absolute error and CPU time (seconds) of Example 3-Eq. (67)
Table 4 Comparison between different methods to solve Eq. (67)
Fig. 2
figure 2

Figures of the approximate solution \(u_{13}(t)\) and the absolute error \(|u-u_{13}|\) for Example  3–Eq. (67)

Example 4

Consider the fourth-order Emden-Flower type equation [32, 36, 37]

$$\begin{aligned} u^{(4)}(t)+\dfrac{3}{t}u^{(3)}(t)+u^{2}(t)-u^{3}(t)=\dfrac{18}{t}+(30+11\, t+t^2)e^{t}+t^{4}e^{2t}-t^{6}e^{3t},\,\,\,\,\,0\le t\le 1, \end{aligned}$$
(68)

subject to the initial-boundary conditions

$$\begin{aligned} u(0)=u^{(1)}(0)=0,\,u(1)=e,\,u^{(1)}(1)=3e, \end{aligned}$$

with the exact solution: \(u(t)=t^{2}e^{t}.\)

Table 6 compares our technique to the methods proposed in [32, 36, 37], while Table 5 shows the highest errors that arise from using our method (SC3GCOMMB2). Further, Fig. 3a shows the computational error norm with \(N = 11\), and Fig. 3b shows the approximate and precise solutions.

Table 5 Maximum absolute error and CPU time (seconds) of Example 4
Table 6 Comparison between different methods of Example 4
Fig. 3
figure 3

Figures of the approximate solution \(u_{11}(t)\) and the absolute error \(|u-u_{11}|\) for Example 4

Example 5

Consider the non-linear Emden-Flower type equation

$$\begin{aligned} u^{(7)}(t)+\dfrac{4}{t}u^{(6)}(t)-(u^{(5)}(t))^{2}-(u(t))^{3}\,= f(t),\,\,\,0\le t\le 1, \end{aligned}$$
(69)

subject to the initial-boundary conditions

$$\begin{aligned} u^{(i)}(0)=0,\,\,i=0,1,2,3,\,\,u^{(i)}(1)=0,\,\,i=0,1,2, \end{aligned}$$

where f(t) is selected to meet the exact solution \(u(t)=t^4(1-t)^3\,e^{t}.\)

The greatest errors that resulted from using our SC3GCOMMB3 technique are shown in Table 7. Figure 4a shows the computational error norm with \(N = 14\), and Fig. 4b shows the approximate and precise solutions.

Table 7 Maximum absolute error and CPU time (seconds) of Example 5
Fig. 4
figure 4

Figures of the approximate solution \(u_{14}(t)\) and the absolute error \(|u-u_{14}|\) for Example 5

Example 6

Consider the non-linear Lane-Emden equation

$$\begin{aligned} u^{(2)}(t)+\dfrac{2}{t}u^{(1)}(t)-e^{-u(t)}=0,\,\,\,0\le t\le 1, \end{aligned}$$
(70)

subject to the initial-boundary conditions

$$\begin{aligned} u(0)=0,\,\,u(1)=0, \end{aligned}$$

where the explicit exact solution is not available, so the following error norm is used to check the accuracy in this case:

$$\begin{aligned} \Vert E_{N}\Vert _{\infty }=\underset{t\in [0,1]}{max}|t\,D^{2}u_{N}(t)+2\,D u_{N}(t)-t\,e^{-u_{N}(t)}|. \end{aligned}$$
(71)

Utilizing SC3GCOMMB2 and using \(N=1,3,6,9,12,14,\) give the numerical results shown in Table 8.

Table 8 Maximum absolute error of Example 6

7 Results and Discussions

New spectral solutions of the singular high-order type equations were implemented and analyzed. Three types of high-order BVPs were solved by applying a unified spectral collocation algorithm. Three modified CPs were utilized as basis functions. The core of the unified algorithm is the establishment of three operational matrices of the modified CPs. The convergence analysis of the three modified Chebyshev expansions was investigated. Several singular BVPs of different orders and various BCs were solved via the proposed technique. The numerical results show the agreement of the resulting approximate solution with the exact ones. As an expected future work, we aim to introduce other modified OPs to treat other types of DEs. The presented method can be developed to solve partial differential equations, such as diffusion-wave and heat equations. The computed operational matrices can also be generalized for solving fractional differential equations.