This example shows how to solve Emden's equation, which is a boundary value problem with a singular term that arises in modeling a spherical body of gas.

After reducing the PDE of the model using symmetry, the equation becomes a second-order ODE defined on the interval $\left[0,1\right]$,

$${y}^{\prime \prime}+\frac{2}{x}{y}^{\prime}+{y}^{5}=0$$.

At $\mathit{x}=0$, the $\left(2/\mathit{x}\right)$ term is singular, but symmetry implies the boundary condition ${\mathit{y}}^{\prime}\left(0\right)=0$. With this boundary condition, the term $\left(2/\mathit{x}\right){\mathit{y}}^{\prime}$ is well defined as $\mathit{x}\to 0$. For the boundary condition $\mathit{y}\left(1\right)=\sqrt{3}/2$, the BVP has an analytical solution that you can compare to a numeric solution calculated in MATLAB®,

$$y\left(x\right)={\left[\sqrt{\frac{1+{x}^{2}}{3}}\right]}^{-1}$$.

The BVP solver `bvp4c`

can solve singular BVPs that have the form

$${y}^{\prime}=S\frac{y}{x}+f(x,y)$$.

The matrix $\mathit{S}$ must be constant and the boundary conditions at $\mathit{x}=0$ must be consistent with the necessary condition $\mathit{S}\cdot \mathit{y}\left(0\right)=0$. Use the `'SingularTerm'`

option of `bvpset`

to pass the `S`

matrix to the solver.

You can rewrite Emden's equation as a system of first-order equations using ${\mathit{y}}_{1}=\mathit{y}$ and ${\mathit{y}}_{2}={\mathit{y}}^{\prime}$ as

${{\mathit{y}}_{1}}^{\prime}={\mathit{y}}_{2}$,

${{\mathit{y}}_{2}}^{\prime}=-\frac{2}{\mathit{x}}{\mathit{y}}_{2}-{\mathit{y}}_{1}^{5}$.

The boundary conditions become

${\mathit{y}}_{2}\left(0\right)=0$,

${\mathit{y}}_{1}\left(1\right)=\sqrt{3}/2$.

In the required matrix form, the system is

$\left[\begin{array}{c}{{\mathit{y}}_{1}}^{\prime}\text{\hspace{0.17em}}\\ {{\mathit{y}}_{2}}^{\prime}\end{array}\right]=\frac{1}{\mathit{x}}\left[\begin{array}{cc}0& 0\\ 0& -2\end{array}\right]\left[\begin{array}{c}{\mathit{y}}_{1}\\ {\mathit{y}}_{2}\end{array}\right]+\left[\begin{array}{c}{\mathit{y}}_{2}\\ -{\mathit{y}}_{1}^{5}\end{array}\right]$.

In matrix form it is clear that $\mathit{S}=\left[\begin{array}{cc}0& 0\\ 0& -2\end{array}\right]$ and $\mathit{f}\left(\mathit{x},\mathit{y}\right)=\left[\begin{array}{c}{\mathit{y}}_{2}\\ -{\mathit{y}}_{1}^{5}\end{array}\right]$.

To solve this system of equations in MATLAB, you need to code the equations, boundary conditions, and options before calling the boundary value problem solver `bvp4c`

. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Create a function to code the equations for $\mathit{f}\left(\mathit{x},\mathit{y}\right)$. This function should have the signature `dydx = emdenode(x,y)`

, where:

`x`

is the independent variable.`y`

is the dependent variable.`dydx(1)`

gives the equation for ${{\mathit{y}}_{1}}^{\prime}$, and`dydx(2)`

the equation for ${{\mathit{y}}_{2}}^{\prime}$.

These inputs are automatically passed to the function by the solver, but the variable names determine how you code the equations. In this case:

function dydx = emdenode(x,y) dydx = [y(2) -y(1)^5]; end

The term that contains `S`

is passed to the solver using options, so that term is not included in the function.

Now, write a function that returns the residual value of the boundary conditions at the boundary points. This function should have the signature `res = emdenbc(ya,yb)`

, where:

`ya`

is the value of the boundary condition at the beginning of the interval.`yb`

is the value of the boundary condition at the end of the interval.

For this problem, one of the boundary conditions is for ${\mathit{y}}_{1}$, and the other is for ${\mathit{y}}_{2}$. To calculate the residual values, you need to put the boundary conditions into the form $\mathit{g}\left(\mathit{x},\mathit{y}\right)=0$.

In this form the boundary conditions are

${\mathit{y}}_{2}\left(0\right)=0$,

${\mathit{y}}_{1}\left(1\right)-\sqrt{3}/2=0$.

The corresponding function is then

function res = emdenbc(ya,yb) res = [ya(2) yb(1) - sqrt(3)/2]; end

Lastly, create an initial guess of the solution. For this problem, use a constant initial guess that satisfies the boundary conditions, and a simple mesh of five points between 0 and 1. Using many mesh points is unnecessary since the BVP solver adapts these points during the solution process.

${\mathit{y}}_{1}=\sqrt{3}/2$,

${\mathit{y}}_{2}=0$.

guess = [sqrt(3)/2; 0]; xmesh = linspace(0,1,5); solinit = bvpinit(xmesh, guess);

Create a matrix for `S`

and pass it to `bvpset`

as the value of the `'SingularTerm'`

option. Finally, call `bvp4c`

with the ODE function, boundary condition function, initial guess, and option structure.

```
S = [0 0; 0 -2];
options = bvpset('SingularTerm',S);
sol = bvp4c(@emdenode, @emdenbc, solinit, options);
```

Calculate the value of the analytical solution in $\left[0,1\right]$.

x = linspace(0,1); truy = 1 ./ sqrt(1 + (x.^2)/3);

Plot the analytical solution and the solution calculated by `bvp4c`

for comparison.

plot(x,truy,sol.x,sol.y(1,:),'ro'); title('Emden Problem -- BVP with Singular Term.') legend('Analytical','Computed'); xlabel('x'); ylabel('solution y');

Listed here are the local helper functions that the BVP solver `bvp4c`

calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function dydx = emdenode(x,y) % equation being solved dydx = [y(2) -y(1)^5]; end %------------------------------------------- function res = emdenbc(ya,yb) % boundary conditions res = [ya(2) yb(1) - sqrt(3)/2]; end %-------------------------------------------

`bvp4c`

| `bvp5c`

| `bvpinit`

| `bvpset`