## 1 Introduction

We consider the following system of singular boundary value problems (SBVPs)

(1.1) |

where are real constants. Here, , , and with .

We next consider the following system of Lane-Emden equations [2, 3, 4, 5, 6, 7, 1] a particular case of (1.1) with and as

(1.2) |

In recent years, the SBVPs have been studied extensively in [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30] and references therein. But as far as we know, we find only the following results on system of Lane-Emden equations. Recently, in [31, 5, 32] authors studied (1.2) with boundary conditions and shape factors and that relates the concentration of the carbon substrate and the concentration of oxygen, the Michaelis-Menten functions

where and are the respective Michaelis-Menten nonlinear operators.

In [3, 4, 33], authors considered the system of Lane–Emden equations (1.2) with boundary conditions and shape factors and occurs in catalytic diffusion reactions with the parameters , are actual chemical reactions.

In [4, 32], the ADM was applied to obtain a convergent analytic approximate solution of (1.2) with . Later, in [34], the variational iteration method (VIM) was applied to obtain approximations to solutions of (1.2) for shape factors . Most recently, a numerical procedure based on sinc-collocation method was developed in [35] to obtain the solution of (1.2). In [36] authors used the reproducing kernel Hilbert space method for solving to obtain the solution of (1.2). In [1], the ADM with Green’s function used to find numerical approximation of the solutions of (1.2).

The homotopy analysis method was developed and improved by S. Liao [37, 38, 39, 40] for solving a broad class of functional equations. Various modifications of the homotopy analysis method have also been elaborated, for example, the optimal homotopy asymptotic method was proposed by Marinca and Herisanu [41, 42], the optimal homotopy analysis method was introduced in [43, 44]. Recently, the homotopy analysis method with Green’s function was applied to solve singular boundary value problems in [45, 46, 47, 48, 49].

In this work, we present the homotopy analysis method combined with the Green’s function strategy for obtaining approximate solutions of coupled singular boundary value problems (1.1) and (1.2). Unlike standard HAM, our approach avoids solving the transcendental equations for the undetermined coefficients. Unlike the ADM and VIM, our proposed technique contains the convergence-control parameters which gives fast convergence of the series solution. Numerical results reveal that the present process provides better results as compared to some existing iterative methods [4, 5, 1].

## 2 Integral form of system of Lane-Emden-Fowler types equations

By following [20], we transform the following system of Lane-Emden-Fowler types equations with Neumann-Robin boundary conditions

(2.1) |

into the equivalent system of integral equations

(2.2) |

where and are given by

(2.3) |

with ,

(2.4) |

with

## 3 Solution method

We consider the system of integral operator equation (2.2) as

(3.1) |

Note that the above system of integral equations reduce to (2.6) with and .

In order to apply the HAM with Green’s function [45, 46, 47, 48, 49], we construct the zeroth-order deformation equations for (3.1) as

(3.2) |

where is an embedding parameter, are initial approximations, are convergence control parameters and () are unknown functions.

If and the homotopy equations (3.2) vary from . Now as varies from to , the solutions of (3.1) will vary from the initial guesses to the exact solutions of equations (3.1).

Expanding and as a Taylor series with respect to yields

(3.3) |

where

(3.4) |

The series (3.3) converge for if are chosen properly, and they take the form

(3.5) |

which will be the solutions of (3.1).

For further analysis, we define the vectors as

(3.6) |

Differentiating (3.2) -times with respect to , dividing them by , setting we obtain the th-order deformation equations as

(3.7) |

where and

On simplification we get

(3.8) |

where and are given by

(3.9) |

Making use of (3.8), the th-order deformation equations (3.7) are simplified as

(3.10) |

By taking the initial approximations and , the solutions components are computed as

(3.11) |

The th-order approximations to solutions are obtained as

(3.12) |

The unknown parameters have a significant influence on the convergence of the approximate solution. We next find the optimal values of by solving

(3.13) |

where and are given by

(3.14) |

Then these optimal values of and will be substituting in (3.12) to obtain the optimal HAM-approximate solution

(3.15) |

## 4 Convergence analysis

Let be a Banach space with norm

(4.1) |

where, and . We next discuss the convergence and error analysis of the proposed method. To do so, we first introduce the vector notations

and

###### Theorem 4.1.

Suppose that the nonlinear function satisfy Lipschitz condition: where and are Lipschitz constants. If the parameter is chosen such that there exists constant , then the series defined by (3.15) is convergent in Banach space .