Singular solutions of a system of nonlinear 2nd order ODEs I'm faced with the following nonlinear 2nd order system of ODEs:
$$
\phi''(r)+\frac{4r^3-1}{r^4-r}\phi'(r)+\frac{r^2 h(r)^2+2r(r^3-1)}{(r^3-1)^2}\phi(r)=0, \\
h''(r)+\frac{2}{r}h'(r)-\frac{2r\phi(r)^2}{r^3-1}h(r)=0.
$$
Here, $\phi$ and $h$ are real-valued functions of $r\in[1,\infty)$ and $h(1)=0$ and $\phi(1)$ is finite and nonzero.
I've been able to find numerical solutions of this system. I did so by assuming $h$ and $\phi$ can be written as a power series around $1$ and by imposing initial conditions at $r=1.00001$ and integrating to larger $r$.
However, I am wondering now if the solutions I found are unique, because the system is singular at the point $r=1$. 
So my question is: for this specific system, is it possible that there exist singular solutions (by which I mean that they are not uniquely fixed by the initial value problem at $r=1.00001$)? Or are the solutions I can find uniquely fixed by the initial values because I impose them at $r=1.00001$ rather than at $r=1$? Note: I am only looking for solutions which in fact do behave analytically at $r=1$!
Any reference to existence theorems relevant for systems like the one above is also appreciated!
 A: I have made a Matlab implementation of your system (script below) to get approximate solutions.
a) If $h(1)=0$ (your hypothesis), function $h$ is found by the solver to be identically zero. This is one of the singular solution you are looking for... even it is rather trivial, isn't it ? Fonctions $\phi$ have are decreasing steadily, in a smooth manner.
I have thus attempted to solve the first equation by replacing $h(r)$ by $0$.
I have turned to Mathematica for solving it;                                                                    
DSolve[p''[r]+(4r^3-1)p'[r]/(r(r^3 - 1))+2 r p[r]/(r^3-1) == 0, p[r], r] 
gives two exact independent solutions : 
$p_1[r]= \ $ Hypergeometric2F1[1/3, 2/3, 1, r^3] and 
$p_2[r]= \ $ MeijerG[{{}, {1/3, 2/3}}, {{0, 0}, {}}, r^3]}} 
(I never met any Meijer before !). 
The general solution is thus a linear combination of these two.
BUT a major drawback is that this solution is defined only for $0<r<1$... Should there exist an analytical continuation through the complex plane ?

b) If one relaxes the constraint on $h(0)$ by letting $h(0)=0.001$ for example, the solutions for $h$ are increasing functions ; the curves of the $phi$ functions are highly oscillatory (see figures which are but a number of very different cases I have tried). 

Matlab program:
function soldiffequ;
clear all;close all;
rspan = 1.00001:0.001:15; % domain for r values 
ode = @(r,y) sc(r,y);
for k=2:5
y0 = [2^k;0.0001;0.;0.];% initial values for phi,h,phi',h', resp.
[r,Y] = ode45(ode, rspan, y0);%Runge Kutta solver
functions phi and h as functions of r
figure(1);hold on;plot(r,Y(:,1),’r’); % phi(r)
figure(2);hold on;plot(r,Y(:,2),’k’); % h(r)
end;
function yp = sc(r,y)
phi=y(1);h=y(2);Phi=y(3);H=y(4);
Dphi=Phi;% a name for the der. of phi
Dh=H;%a name for the der. of h
DPhi=-((4*r.^3-1)./(r.^4-r))Phi-phi.((r.^2.*h.^2+2*r*(r.^3-1))./((r.^3-1).^2));
DH=-(2./r).H+h.((2*r.*(phi.^2))./(r.^3-1));
yp = [Dphi;Dh;DPhi;DH];

Fig. 1a : case h(1)=0 : in red, functions phi for initial values 4,8,16 ,32.
Functions h remain identically 0… :
Fig. 2a and 2b : Cases  where h(1) is non zero (though very small : 0.0001) showing that functions phi may have a very oscillatory behaviour, according to their initial values  (here 4,8,16,32),; 
Functions h are no longer zero : they are even rapidly increasing.



