Yeah, sure, I forget the x'(0) initial cond. However I don't see how it is possible to write the system with onlye 2 state variables, but I think I figured out how to do it with 3 state variables. Thanks

Yeah, you are right, my fault, but the notations are x1 = x, x2 = x' = x1'.

The ODE follows from the Lagrangian

L=x^(-2/a)*x'^2-b*x^(-2/a).

Therefore you have the Hamiltonian constrain

h=x^(-2/a)* (x'^2+b) where h=x(0)^(-2/a)*(x'(0)^2+b)

Hence x'=sqrt(h*x^(2/a)-b).

You can solve this problem analytically.

If you let dx/dt = y, then dy/dt = (y^2 + b)/(a*x)

then x(0) = x0 and x'(0) = y0

You can write: (dy/dt) / (dx/dt) = dy/dx = (y^2 + b)/(a*x*y)

This is separable: y/(y^2 + b) * dy = dx/(a*x)

Integrating: 1/2*ln( (y^2 + b) / (y0^2 + b) ) = 1/a*ln(x/x0)

Removing the "ln": ( (y^2 + b) / (y0^2 + b) )^(1/2) = (x/x0)^(1/a)

Solve for "x": x = x0*( (y^2 + b) / (y0^2 + b) )^(a/2)

and substitute into the original equation:

dy/dt = (y^2+b)^(1-a/2) * (y0^2 + b) )^(a/2) / (a*x0)

= c*(y^2+b)^(1-a/2)

Where c = (y0^2 + b) )^(a/2) / (a*x0)

The ODE is again separable:

[(y^2+b)^(a/2-1)] * dy = [c] * dt

The LHS is integrable, but the form of the function depends on the value of "a".

You will end up with f(y) = f(y0) + ct and x = (y^2 + b) / (a*c)

Example 1: a = 3, b = 0

y = sqrt(2*c*t + y0^2)

x = x0*( y/y0 )^3

Example 2: a = 2

y = c*t + y0

x = x0* (y^2 + b) / (y0^2 + b)

For random values of "a", it may get tricky.