**4. Method of Solution**

Approximating the unknown function in terms of orthogonal polynomials has been practiced in several papers in recent years [18,21,22,32,33] for different types of problems. Here, for solving the problem in Equation (1), we approximate

$$D^{\mathfrak{a}}y(\mathfrak{x}) = \mathbb{C}^{T}\Phi\_{\mathfrak{n}}(\mathfrak{x}).\tag{19}$$

We are approximating the derivative first because we want to use the initial condition. Taking the integral of order *α* on both sides of Equation (19), we get

$$y(x) = \mathbb{C}^T I^a \Phi\_n(x) + y(0). \tag{20}$$

Using the operational matrix of integration, Equation (20) can be written as

$$\Phi(\mathbf{x}) \cong \mathbb{C}^{T} I^{(a)} \Phi\_{\mathfrak{n}}(\mathbf{x}) + A^{T} \Phi\_{\mathfrak{n}}(\mathbf{x}) \tag{21}$$

where *y*(0) = *a* ∼= *AT*Φ*n*(*x*) and *I*(*α*) is the operational matrix of integration of order *α*.

Using Equation (19), we can write

$$\mathcal{I}I^{1-a}y(\mathbf{x}) = I\mathcal{D}^a y(\mathbf{x}) = \mathcal{C}^T I\Phi\_\mathbb{R}(\mathbf{x}) \cong \mathcal{C}^T I^{(1)}\Phi\_\mathbb{R}(\mathbf{x}).\tag{22}$$

Using Equations (19) and (22) in Equation (1), we obtain

$$J\left(\mathbf{c}\_{0},\mathbf{c}\_{1},\ldots,\mathbf{c}\_{n}\right) = \int\_{0}^{1} \left(\mathbf{g}(\mathbf{x})\mathbf{C}^{T}\boldsymbol{\Phi}\_{\text{n}}(\mathbf{x}) + \mathbf{g}'(\mathbf{x})\mathbf{C}^{T}I\boldsymbol{\Phi}\_{\text{n}}(\mathbf{x}) + h'(\mathbf{x})\right)^{2}d\mathbf{x}.\tag{23}$$

Equation (23) can then be written as

$$J\left(\mathbf{c}\_{0},\mathbf{c}\_{1},\ldots,\mathbf{c}\_{n}\right) = \int\_{0}^{1} \left(\mathbb{C}^{T}\mathbf{g}(\mathbf{x})\Phi\_{\mathbf{n}}(\mathbf{x}) + \mathbb{C}^{T}I^{(1)}\mathbf{g}'(\mathbf{x})\Phi\_{\mathbf{n}}(\mathbf{x}) + h'(\mathbf{x})\right)^{2}d\mathbf{x}.\tag{24}$$

We further take the following approximations:

$$\mathcal{S}\_{\mathcal{S}}(\mathbf{x})\Phi\_i(\mathbf{x}) \cong E\_1^{i,T}\Phi\_n(\mathbf{x})\tag{25}$$

$$\log'(\mathbf{x})\Phi\_i(\mathbf{x}) \cong E\_2^{i\_T T} \Phi\_n(\mathbf{x})\tag{26}$$

$$h'(\mathbf{x}) \cong E\_3^T \Phi\_n(\mathbf{x}) \tag{27}$$

where *Ei*,*<sup>T</sup>* <sup>1</sup> = [*e<sup>i</sup>* 1,0, *<sup>e</sup><sup>i</sup>* 1,1, ...,*e<sup>i</sup>* 1,*n*], *<sup>E</sup>i*,*<sup>T</sup>* <sup>2</sup> = [*e<sup>i</sup>* 2,0, *<sup>e</sup><sup>i</sup>* 2,1, ...,*e<sup>i</sup>* 2,*n*], *<sup>E</sup><sup>T</sup>* <sup>3</sup> = [*e*3,0, *e*3,1, ...,*e*3,*n*], and *ei* 1,*<sup>j</sup>* <sup>=</sup> *g*(*x*)Φ*i*(*x*), <sup>Ψ</sup>*j*(*x*), *<sup>e</sup><sup>i</sup>* 2,*<sup>j</sup>* = *g* (*x*)Φ*i*(*x*), Ψ*j*(*x*), *e*3,*<sup>j</sup>* = *h* (*x*), Ψ*j*(*x*), 0 ≤ *i*, *j* ≤ *n*, and −, − is the usual inner product space.

Using Equations (25) and (26) we can write

$$\Phi\_{\mathcal{S}}(\mathbf{x})\Phi\_{\mathfrak{n}}(\mathbf{x}) \cong E\_1^T \Phi\_{\mathfrak{n}}(\mathbf{x}) \tag{28}$$

$$g'(\mathbf{x})\Phi\_{\mathbf{n}}(\mathbf{x}) \cong E\_2^T \Phi\_{\mathbf{n}}(\mathbf{x})\tag{29}$$

where

$$E\_1^T = \left(E\_1^{i,T}\right)\_{0 \le i \le n} \text{ and } E\_2^T = \left(E\_2^{i,T}\right)\_{0 \le i \le n}. \tag{30}$$

From Equations (24) and (27)–(29), we get

$$J\left(\mathbf{c}\_{0},\mathbf{c}\_{1},\ldots,\mathbf{c}\_{n}\right) = \int\_{0}^{1} \left(\mathbf{C}^{T}E\_{1}^{T}\Phi\_{n}(\mathbf{x}) + \mathbf{C}^{T}I^{(1)}E\_{2}^{T}\Phi\_{n}(\mathbf{x}) + E\_{3}^{T}\Phi\_{n}(\mathbf{x})\right)^{2}d\mathbf{x}.\tag{31}$$

Let

$$E^T = \mathbb{C}^T \left( E\_1^T + I^{(1)} E\_2^T \right) + E\_3^T. \tag{32}$$

From Equations (31) and (32), we get

$$\begin{array}{rcl} \int \left( \mathbf{c}\_{0}, \mathbf{c}\_{1}, \dots, \mathbf{c}\_{n} \right) & = \int\_{0}^{1} \left( E^{T} \Phi\_{n}(\mathbf{x}) \right)^{2} d\mathbf{x} \\ & = \int\_{0}^{1} E^{T} \Phi\_{n}(\mathbf{x}) \Phi\_{n}(\mathbf{x})^{T} E \, \mathbf{d} \mathbf{x}, \\ & = E^{T} P E \end{array} \tag{33}$$

where *P* is a square matrix given by *P* = - 1 <sup>0</sup> Φ*n*(*x*)Φ*n*(*x*) *<sup>T</sup> dx*.

Using Equation (22), the boundary condition can be written as

$$I^{1-a}y(1) \cong \mathbb{C}^T I^{(1)} \Phi\_n(1) = \epsilon. \tag{34}$$

Using the Lagrange multiplier method [18,20–22,32,33], the necessary extremal condition for the functional in Equation (33) becomes

$$
\frac{\partial f}{\partial \mathcal{L}\_0} = 0, \ \frac{\partial f}{\partial \mathcal{L}\_1} = 0, \dots, \ \frac{\partial f}{\partial \mathcal{L}\_{n-1}} = 0. \tag{35}
$$

From Equations (34) and (35), we get a set of *n* + 1 equations. Solving these *n* + 1 equations, we get unknown parameters *c*0, *c*1, ... , *cn*. Using these unknown parameters in Equation (21), we get the unknown function's extreme values of the non-linear fractional functional.
