Question: Error : unable to convert to explicit first-order system

Hello,

I am trying to solve the following differential coupled system :

DN:=diff(N(x),x)=a/4*(x+8*(sin(F(x)))^2/x)*N(x)*diff(F(x),x)^2;

Du:=diff(u(x),x)=a/8*((x^2+8*(sin(F(x)))^2)*(1-u(x)/x)*diff(F(x),x)^2+2*(sin(F(x)))^2+4*(sin(F(x)))^4/x^2);

DDF:=diff(F(x),x$2)=1/(N(x)*(1-u(x)/x)*(x^2+8*(sin(F(x)))^2))*(-(x^2+8*(sin(F(x)))^2)*diff(N(x),x)*(1-u(x)/x)*diff(F(x),x)+(1+4*(sin(F(x)))^2/x^2+4*(1-u(x)/x)*diff(F(x),x)^2)*N(x)*sin(2*F(x))-2*(x+4*sin(2*F(x))*diff(F(x),x))*N(x)*(1-u(x)/x)*diff(F(x),x)-2*(1+8*(sin(F(x)))^2/x^2)*(u(x)-diff(u(x),x)*x)*N(x)*diff(F(x),x));

I developped a series solution around 0 in order to have some boundary conditions just next to the singularity at x=0 :

f0(x):=Pi+b1*x;

df0(x):=diff(f0(x),x);

u0(x):=a/8*b1^2*(1+4*b1^2)*x^3;

du0(x):=diff(u0(x),x);

n0(x):=b2+a/8*b1^2*b2(1+8*b1^2)*x^2;

dn0(x):=diff(n0(x),x);

I have some 11 boundary conditions at x1 (small) and x2 (big):

x1:=0.001;

x2:=5000;

f1:=eval(f0(x),x=x1);

u1:=eval(u0(x),x=x1);

n1:=eval(n0(x),x=x1);

f_1:=eval(df0(x),x=x1);

u_1:=eval(du0(x),x=x1);

n_1:=eval(dn0(x),x=x1);

f2:=0;

f_2:=0;

n2:=1;

n_2:=0;

u_2:=0;

I put all of these in a system and tried to solve using dsolve with the middefer bvp method (I tried others too) :

dsys:={DN,Du,DDF,F(x1)=f1,u(x1)=u1,N(x1)=n1,diff(F(x),x)(x1)=f_1,diff(u(x),x)(x1)=u_1,diff(N(x),x)(x1)=n_1,F(x2)=f2,diff(F(x),x)(x2)=f_2,N(x2)=1,diff(u(x),x)(x2)=0,diff(N(x),x)(x2)=0};

dsoln:=dsolve(dsys,numeric,method=bvp[middefer]);

Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system

I also had to give some values to b1 and b2, since there was an error saying I had too many independent variables. Those values, to my knowledge, are good estimations :

b1:=-1;

b2:=1;

What to do with this? How can I get dsolve to solve this system? I have been working on this for a while now, and also checked some possible solutions in the other questions. I noticed that having terms in (dF/dx)^2 could be a problem, so i rewrote my equations :

reweq1:=diff(F(x),x)=sqrt(4*diff(N(x),x)/(a*N(x)*(x+8*(sin(F(x)))^2/x)));

reweq2:=diff(F(x),x)=sqrt((8*diff(u(x),x)/a-2*(sin(F(x)))^2-4*(sin(F(x)))^4/x^2)/((1-u(x)/x)*(x^2+8*(sin(F(x)))^2)));

reweq3:=diff(F(x),x)=sqrt((-1-4*(sin(F(x)))^2/x^2+(diff(F(x),x$2)*N(x)*C*(x^2+8*(sin(F(x)))^2)+(x^2+8*(sin(F(x)))^2)*diff(N(x),x)*C*diff(F(x),x)+2*(x+4*sin(2*F(x))*diff(F(x),x))*N(x)*C*diff(F(x),x)+2*(1+8*(sin(F(x)))^2/x^2)*(u(x)-diff(u(x),x)*x)*N(x)*diff(F(x),x))/(N(x)*sin(2*F(x))))/4/C);

Same error, doesn't work.

Thanks in advance for any help!

PS:For the curious among you, those equations come from the Skyrme-Einstein model for black holes :)