Solving the ODE

(* Mathematica does not know how to solve the equation *)

DSolve[y''[x]==y[x]^(3/2)/Sqrt[x],y,x]

(* We can try to solve it numerically by starting from eps to avoid the singularity at 0 *)
eps=10^-3
eqn=NDSolve[{y''[x]==x^-0.5 y[x]^1.5,y[eps]==1,y'[eps]==-2},y,{x,eps,10}]
{{xmin,xmax}}=(y/.eqn[[1]])["Domain"]
Plot[y[x]/.eqn,{x,xmin,xmax}]

(* Naturally if we start from x_start = eps, then y(x_start) is not exactly 1. A first order correction could improve the situation? *)
ϵ=10^-6;
c=-1.588071 (* e.g., see https://en.wikipedia.org/wiki/Thomas%E2%80%93Fermi_equation *)
y0=c ϵ+1
eqn=NDSolve[{y''[x]==x^-0.5 y[x]^1.5,y[ϵ]==y0,y'[ϵ]==c},y,{x,ϵ,5}]
{{xmin,xmax}}=(y/.eqn[[1]])["Domain"]
des0=Plot[y[x]/.eqn,{x,xmin,xmax}]

(* Naturally if we start from x_start = eps, then y(x_start) is not exactly 1. A first order correction could improve the situation? *)
ϵ=10^-4;
c=-1.56841
y0=c ϵ+1
eqn=NDSolve[{y''[x]==x^-0.5 y[x]^1.5,y[ϵ]==y0,y'[ϵ]==c},y,{x,ϵ,8}]
{{xmin,xmax}}=(y/.eqn[[1]])["Domain"]
des1=Show[Plot[y[x]/.eqn,{x,xmin,xmax},PlotStyle->Red],des0]

Phase portrait

(* phase portrait for the associated differential equation *)
Show[StreamPlot[{1,(4q+p^1.5)/(3p+q)},{p,0,200},{q,0,-500}],
Plot[+p^1.5/-4,{p,0,200},PlotStyle->Red],
Plot[-3p ,{p,0,200},PlotStyle->Black]]

Asymptotic behaviour

(* the asymptotic behaviour is determined as one of the singular points *)

Solve[{4q+p^1.5==0,3p+q==0},{p,q}]

(* the asymptotic behaviour can also be determined by substituting y = A x^-3 into the differential equation *)

eq=ReleaseHold[HoldForm[D[y[x],x,x]==y[x]^(3/2)/Sqrt[x]]/.y[x]->A x^-3]

FullSimplify[eq,Assumptions->{A>0,x>0}]

Solving the associated ODE

(* we apply the L'Hôpital rule for finding the tangent at singular point P, and we select the one that corresponds to the sepatrix *)
eq=Solve[m==(4m+18)/(3+m),m]
mm=m/.eq[[1]]

(* we now solve the associated differential equation by starting from epsilon from P to avoid the singularity *)
ϵ=10^-7;
p0=144-ϵ;
eqn=NDSolve[{
(4q[p]+p^1.5)/(3p+q[p])==q'[p],
q[p0]==-432-mm ϵ},q,{p,p0,ϵ},AccuracyGoal->∞]

des=Plot[q[p]/.eqn,{p,0,p0},PlotStyle->Purple]

(* we plot the separatrix in the phase plane *)

Show[
Show[StreamPlot[{1,(4q+p^1.5)/(3p+q)},{p,0,200},{q,0,-500}],
Plot[+p^1.5/-4,{p,0,200},PlotStyle->Red],
Plot[-3p ,{p,0,200},PlotStyle->Black]],
des]

(* we now solve the associated differential equation by starting from epsilon from P to avoid the singularity *)
ϵ=10^-7;
p0=144-ϵ;
eqnT=NDSolve[{
(4q[p]+p^1.5)/(3p+q[p])==q'[p],
q[p0]==-432-mm ϵ},q,{p,p0,0}  ,AccuracyGoal->∞,MaxSteps->Infinity]

(* Note that q/p^(4/3) gives an approximation of y'(0) *)
val=(q[#]/#^(4/3)/.eqnT[[1]])&/@{10^-1,10^-3,10^-6,10^-9,10^-12,10^-15}
(* and this correspond to a value of x = q/p/y' *)
((q[#]/# /.eqnT[[1]])&/@{10^-1,10^-3,10^-6,10^-9,10^-12,10^-15})/val

Integration of the separatrix

fnc[z_?NumericQ]:=q[z]/.eqnT[[1]]
Plot[fnc[p],{p,0,144}]

(* if we integrate from P... hum this does gives us y = 144 / x^3 *)
xlim=100;
eqnN=NDSolve[{y'[x]==fnc[x^3 y[x]]/x^4,y[xlim]==144/xlim^3},y,{x,xlim,0.1},AccuracyGoal->∞]
 y[2]/.eqnN 
144/2^3

(* if we integrate from 0... sounds better *)
xlim=10^-4;
eqnN2=NDSolve[{y'[x]==fnc[x^3 y[x]]/x^4,y[xlim]==1},y,{x,xlim,20},AccuracyGoal->∞]

 y[2]/.eqnN2
144/2^3

(* comparison of the solutions *)
 Show[ Plot[ 144/x^3,{x,6,20},PlotRange->All,PlotStyle->Dashed], Plot[ y[x]/.eqnN2,{x,1,20},PlotRange->All,PlotStyle->Black],des1]

Show[
Show[StreamPlot[{1,(4q+p^1.5)/(3p+q)},{p,0,5},{q,0,-5}],
Plot[+p^1.5/-4,{p,0,200},PlotStyle->Red],
Plot[-3p ,{p,0,200},PlotStyle->Black]],
des]