Solving a second-order nonlinear ODE with a singularity on x=0 I'm doing some reasearch on electromagnetic nanostructures and I have to solve this differential equation (the exact values of the constants don't matter, I just want all the possible solutions of y(x) given some values to these constants).
$$
\frac{d^2 y}{dx^2}=-\frac{1}{x}\frac{dy}{dx}+\frac{\sin(2y)}{2} (\frac{1}{x^2}+\frac{K}{A})-\frac{D}{A}\frac{\sin(y)}{x}+\frac{\mu HM}{2A}\sin(y)
$$
from x=0 till x=R, with the boundary conditions
$$
y(0)=0,\ \frac{dy}{dx}(R)=\frac{-D}{2A}
$$
I believe you can not find an analitic solution to this equation, so I've been trying to use numerical methods like the shooting method (given the boundary conditions, I found it appropiate).
The thing is that the singularity on x=0 doesn't let me find the solutions. I obtain different results depending on how many steps I take in the method.
I also posted this on Computational Science and Physics StackExchange, but for now I couldn't fix it.
 A: The following maxima script gives you the formulas for the calculation of the coefficients $y_2,\ldots,y_4$ in the ansatz $y=y_1 x + y_2 x^2 + y_3 x^3 + y_4 x^4$. $y_1$ is the free parameter for matching the initial condition at $x=R$.
declare(x,mainvar);

y : y1*x + y2*x^2 + y3*x^3 + y4*x^4 + y5*x^5;
fDiff : x^2 * diff(y,x,2) + x*diff(y,x,1);
fRHS : sin(2*y)/2 + x^2*K/2/A*sin(2*y) - x*D/A*sin(y) + x^2*mu*H*M/(2*A)*sin(y);

eq : factorout(taylor(fRHS,x,0,4) - fDiff,x);

coeff(eq,x,1);
eq1 : A*scsimp(solve(coeff(eq,x,2));
eq2 : 6*A*scsimp(coeff(eq,x,3));
eq3 : 6*A*scsimp(coeff(eq,x,4));

tex(solve([eq1,eq2,eq3],[y2,y3,y4]));

The result is:
\begin{align*}
 {\it y_2}&=-{{D\,{\it y_1}}\over{3\,A}}\\
 {\it y_3}&=-{{4\,A^2\,{\it y_1}^3+\left(A\,\left(-3\,H\,M\,\mu-6\,K\right)-2\,D^
 2\right)\,{\it y_1}}\over{48\,A^2}}\\
 {\it y_4}&={{44\,A^2\,D\,
 {\it y_1}^3+\left(A\,D\,\left(-11\,H\,M\,\mu-22\,K\right)-2\,D^3
 \right)\,{\it y_1}}\over{720\,A^3}}
\end{align*}
You can use the series expansion for $y(x)$ to flee from the singularity and let the numerical solver kick in at some sufficiently large $x$.
You can also play with the order of the Taylor expansion to get more exact solutions or to estimate the quality of the solution for the chosen $x$.
Calculating $\frac{\partial y(x)}{\partial y_1}$ for the compuation of the system matrix of the variational system which you need for the overall-shooting method should be no problem.
The expansion order is an adjustable parameter in the following variant of the above maxima script. This degrades a bit the clarity of the script but it greatly simplifies playing around with the expansion order.
declare(x,mainvar);

order : 4;

y : sum(concat('y,i)*x^i,i,1,order);
fDiff : ratsimp(x^2 * diff(y,x,2) + x*diff(y,x,1));
fRHS : sin(2*y)/2 + x^2*K/2/A*sin(2*y) - x*D/A*sin(y) + x^2*mu*H*M/(2*A)*sin(y);

eq : factorout(taylor(fRHS,x,0,order) - fDiff,x);

for i:1 thru order do concat('eq,i) :: coeff(eq,x,i);

sol : solve(makelist(concat('eq,i),i,2,order),makelist(concat('y,i),i,2,order));

tex(sol);

A: Have you tried linearization? If $y$ is small we can do the following:
Let $\alpha = \frac{K}{A}$, $\beta = \frac{D}{A}$, $\gamma = \frac{\mu HM}{2A}$ and $\delta = -\frac{D}{2A}$ then using $\sin(2y) = 2\sin(y)\cos(y)$ the whole equation can be rewritten as
$$
xy' + x^2 y''  
=
\big((1+\alpha x^2)\underbrace{\cos(y)}_{\approx 1} - \beta x + \gamma x^2 \big)\underbrace{\sin(y)}_{\approx y} 
\simeq
\underbrace{\big( (\alpha + \gamma) x^2 - \beta x +1\big)}_{:=p(x)}y
$$
Then rewriting as a 2 dimensional system of order 1 yields:
$$
\begin{bmatrix}  0 & 1 \\ x & x^2  \end{bmatrix}
\cdot
\frac{d}{dx}
\begin{bmatrix}  y \\ z  \end{bmatrix}
=
\begin{bmatrix}  0 & 1 \\ p(x) & 0 \end{bmatrix}
\cdot
\begin{bmatrix}  y \\ z  \end{bmatrix}
,\qquad 
y(0) = 0, z(R) = \delta
$$
Thus we have rewritten the system as a non-linear DAE $E(x)y' = A(x)y$ with $\det(E) = -\frac{1}{x}$. So here you can  invert $E$:
$$
\frac{d}{dx}
\begin{bmatrix}  y \\ z  \end{bmatrix}
=
\begin{bmatrix}  p(x) & -x \\ 0 & \frac{1}{x} \end{bmatrix}
\cdot
\begin{bmatrix}  y \\ z  \end{bmatrix}
,\qquad 
y(0) = 0, z(R) = \delta
$$
And now you can pass this ODE to your solver, and actually it should even be analytically sovleable.
