# How to reformulate a 1/x^2 singular term to 1/x so that bvp4c can solve it?

I have a ''cosmetic'' problem with a singular term in my Matlab script. I am trying to solve the following system of differential equations: \begin{aligned} y_1' &= y_2,\\ y_2' &= \frac{2y_1}{x^2}+k(y_1^3 - y_1) - \frac{2y_2}{x}. \end{aligned}

Where x is the independent variable of y and k is a constant. This system is to be solved with the boundary conditions: \begin{aligned} y_1(0)&=0,\\ y_1(10) &= 1. \end{aligned}

I have a working script for this with Matlab's bvp4c function. However, the solution for $y_2$ appears to go off to infinity near zero. I've read from Matlab's documentation that bvp4c can solve problems with singular boundary values if they are of the form: $$y' = S \cdot \frac{y}{x} + F(x,y,p).$$ But as you can see there is a $\frac{1}{x^2}$-term in addition to a $\frac{1}{x}$-term. I've been trying to modify the form of the differential equation such that I could write it in the form Matlab wants it to be. However, whatever I do I can't seem to be able to modify it to the desired form. I was hoping that maybe someone might have some insight to this problem. Thanks in advance!

• Have you tried a subsitution, e.g. $Y_{1,2}=\frac{y_{1,2}}{x}$? That should transform the equation into the required form. – nluigi Jan 8 '17 at 14:55
• Yes, but I can't seem to get anywhere with it. Could you show how you would continue with that substitution? I end up with a matrix S of the form: $S = \begin{bmatrix} 0 & 2 \\ 2 & 0 \end{bmatrix}$ Then for the ODE-function I have : $Y_1' = Y_2$ $Y_2' = k(Y_1^3 - Y_1) - 2Y_2$ With the singular term inserted as instructed in Matlab's documentation. – user22872 Jan 8 '17 at 18:15
• Perhaps it might be easier to use a shooting method? It's easy to get the Taylor series at $x=0$ to specify an initial condition at a small positive $x$, and it avoids the problems with bvp4c. – Kirill Jan 8 '17 at 21:19

Your equation, written as a 2nd-order ODE, is $$y'' = -\frac{2}{x}y' + \frac{2}{x^2}y + ky(y^2-1).$$ Substituting $$y(x) = x w(x),$$ it becomes $$w'' = -\frac{4}{x}w' + k w(1 - x^2 w^2),$$ which matches the form accepted by bvp4c with the matrix $$S = \begin{pmatrix} 0 & 0 \\ -4 & 0 \end{pmatrix}.$$ Other transformations of the form $y=(x+\alpha x^{-2})w$ also get it to this kind of form.
There is, in fact, a well-known theory of changes of variables in linear ODEs (see, for example, this theorem due to Liouville: https://mathoverflow.net/a/187499), and finding a transformation like $y=xw$ doesn't actually require any guesswork. Applying the transformations to the first-order form of the ODE in $(y_1,y_2)$ makes this unnecessarily unwieldy.
• Is it a problem that with this transform the BC become $w(0)=\infty$? – nluigi Jan 9 '17 at 18:07
• @nluigi How do you mean? $y(0)=0$ becomes the condition that $w(0)$ is finite, and $w(0)$ is just the unknown $y'(0)$. – Kirill Jan 9 '17 at 18:24
• it seems to me $w(0)=y(0)/0$ is undefined, hence my question. Why do you say it is finite? – nluigi Jan 9 '17 at 19:00
• @nluigi Expand $y$ in Taylor series at $x=0$, then expand $w$ in Taylor series, then equate coefficients in $y=xw$. Writing $w=y/x$ is the wrong way to go about doing that—it's the expression $y(0)/0$ that is undefined, not the boundary condition itself. And in any case, if $y$ is a solution, then $y(0)=0$ and $y$ is a suitably differentiable function, so $w$ definitely exists and is finite. And if $w$ is a solution that remains finite at $0$, then $y=xw$ is a solution of the equation in $y$. – Kirill Jan 9 '17 at 19:02
• @nluigi In fact, a solution to $0=0x$ is not $x=0/0$, which is meaningless, but $x\in\mathbb{R}$. – Kirill Jan 9 '17 at 19:05