Mathematical overview
You seem to be using NDSolve perfectly well, and this differential equation does not seem to pose any special problems for NDSolve. There seems to be two issues raised, the oscillatory solutions and the square root solutions.
It is not uncommon for differential equations to have an isolated non-oscillatory solution among oscillatory ones. With nonlinear equations it is even less clear what the "normal" situation might be. However, we can make some rough guess at how this ODE works. We might view the differential equation in the form
$$(1)\ \ y''(t) + \left(\frac{w^2}{t^2}-\frac{1}{y(t)^4}\right)\;y(t)=0, \quad{\rm or}\quad
(2)\ \ y''(t) + \frac{w^2}{t^2}\,y(t)=\frac{1}{y(t)^3}\,,$$
as a perturbation of the equidimensional equation
$$y''(t) + \frac{w^2}{t^2}\;y(t)=0\,,$$
whose solution is $y(t)=A \sqrt{t} \sin \left({\sqrt{4w^2-1} \over 2}\;\log(t) - t_0\right)$ if $w>1/2$. This in turn can be related to both $y''+ε\,y=0$ and $y''+y/(4t^2)=0$ (i.e. $w=1/2$), The first has periodic solutions if $ε>0$ and exponential growth (generically) if $ε<0$. The second has $y=a\sqrt{t}$ as a solution. If in (1) the coefficient $w^2t^{-2}-y^{-4}$ is greater than $ε>0$ over a long enough interval, then the solution must eventually decrease; and if it is less than $ε<0$ over a long enough interval, the solution must eventually increase. Further if the coefficient can remain close to $1/(4t^2)$, the solution might remain close to a square root or a logarithm times a square root.
When $y(t)$ is large, then the term $-1/y(t)^4$ becomes negligible and the coefficient will be positive. That leads to $y(t)$ decreasing. As it does, the term $-1/y(t)^4$ becomes dominant and the coefficient becomes negative.. Eventually, $y(t)$ must begin to increase. Thus oscillatory behavior seems likely to be a characteristic of many solutions. The way oscillations might be avoided is if in the transition between $y(t)$ large and $y(t)$ small, the solution does not pass across from one region to another but the coefficient approaches zero asymptotically. (It is clear from NDSolve that this does not happen.)
The ODE
Let us set up the ODE this way:
odelincoeff = (w^2/(t^2) - (y[t])^(-4));
ode = y''[t] + odelincoeff y[t]
Power function solutions
The coefficients allow the possibility that there is a solution that is a power function or a sum of power functions. To look for a power function solution, we have to find a power that works, if any. Substitute an arbitrary power function a t^n into the differential formula. The condition for a such a solution is that the result must vanish.
solcondition = ode /. y -> (a #^n &) // Expand
(*
-a n t^(-2 + n) + a n^2 t^(-2 + n) - t^(-3 n)/a^3 + a t^(-2 + n) w^2
*)
The powers must match up in order to get cancellation. That puts a condition on n that we can use to find n.
powercondition = Simplify[
First @ Solve[Equal @@ DeleteDuplicates @ Cases[solcondition, t^_, Infinity], n, Reals],
t > 0]
(*
{n -> 1/2}
*)
Having found a suitable power, we can solve for the undetermined coefficient a:
solcondition /. powercondition
(*
-(1/(a^3 t^(3/2))) - a/(4 t^(3/2)) + (a w^2)/t^(3/2)
*)
coeffcondition = First@Select[
Solve[0 == solcondition /. powercondition /. t -> 1, a],
Positive[a /. # /. w -> 1] &]
(*
{a -> Sqrt[2]/(-1 + 4 w^2)^(1/4)}
*)
Thus a particular solution is this:
powersol = Function[t,
Evaluate[a t^n /. powercondition /. coeffcondition // ToRadicals]]
(*
Function[t, (Sqrt[2] Sqrt[t])/(-1 + 4 w^2)^(1/4)]
*)
For this solution the terms of coefficient of y[t] match degrees and equal the ratio -y''[t]/y[t], which of course it has to.
odelincoeff /. y -> powersol // Simplify
(*
1/(4 t^2)
*)
Note that the coefficient is postive but approaches zero. So while the solution is concave down, it does not, for a fixed [CurlyEpsilon]>0, remain greater than [CurlyEpsilon]. Therefore, unlike the case in the general discussion above, it is not necessary for it to become decreasing, which of course it does not do.
Nearby solutions
To examine what happens near the square root solution, let us first make a change of variables, the standard substitution t = e^s for an equidimensional equation to get a linear differential form plus a nonlinear term -E^(2 s)/u[s]^3.
odelog = (ode /. y -> (u[Log[#]] &) // Factor // Numerator) / u[s]^3 /.
{Log[t] -> s, t -> E^s} // Expand
(*
-(E^(2 s)/u[s]^3) + w^2 u[s] - u'[s] + u''[s]
*)
We can check that our power solution still works.
odelog /. {u -> (powersol[E^#] &)} // Simplify
(*
0
*)
We can examine solution nearby the square root solution by perturbing the solution slightly by a factor (1 + ε v[s]). Replacing ε -> 0 eliminate quadratic and higher order terms in v[s]. (It turns out that the zero-order terms cancel and the ε in the linear terms is canceled, too, by Simplify.)
odeperturbed = Simplify[
0 == (odelog /. {u -> ((1 + ε v[#]) powersol[E^#] &)} // Together // Numerator),
ε != 0] /. ε -> 0
(*
4 (-1 + 4 w^2) v[s] + 4 (v′′)[s] == 0
*)
DSolve[odeperturbed, v[s], s]
(*
{{v[s] -> E^(s Sqrt[1 - 4 w^2]) C[1] + E^(-s Sqrt[1 - 4 w^2]) C[2]}}
*)
Here we convert the exponential function to trigonometric functions. We see that as a function of s = Log[t] the solution is periodic of period2 π / Sqrt[4 w^2-1], but as a function of t, the oscillations will lengthen ast` increases.
vsol = v ->
Function[s,
Evaluate[Simplify[
v[s] /. First @ DSolve[odeperturbed, v[s], s] // ComplexExpand // ExpToTrig,
w > 1 && t > 0] /.
{C[1] -> ε (Cos[s0] + I Sin[s0])/2, C[2] -> ε (Cos[s0] - I Sin[s0])/2} // Simplify]
]
(*
v -> Function[s, ε Cos[s0 + s Sqrt[-1 + 4 w^2]]]
*)
We can check that v makes a good approximation. The constant and linear terms in ε vanish, so for small ε, it will be a good approximation over some interval.
Series[
Simplify[
odelog /. {u -> ((1 + v[#]) powersol[E^#] &)} /. vsol /. {s0 -> 0},
w > 1/2 && s ∈ Reals],
{ε, 0, 2}]
% /. s -> Log[t]
(*
-((3 (E^(s/2) (-1+4 w^2)^(3/4) Cos[s Sqrt[-1+4 w^2]]^2) ε^2)/Sqrt[2]) + O[ε]^3
-((3 (Sqrt[t] (-1+4 w^2)^(3/4) Cos[Sqrt[-1+4 w^2] Log[t]]^2) ε^2)/Sqrt[2]) + O[ε]^3
*)
Thus the solutions near the square root will oscillate back and forth around the square root graph. Pretty cool.
Initial values for the square root
Once one obtains the correct coefficient $a$, the initial values for starting at $t = t_1$ are simply
$$ y(t_1) = a \sqrt{\vphantom{f} t_1}, \quad y'(t_1) =
{a \over 2 \sqrt{\vphantom{f} t_1}}\,.$$
This works for any $t_1>0$, including t0 in the question.
Numerical solution
The function sol gives the solution to the specified initial value problem. The functions rtIvals and rtSol
Clear[sol, rtIvals, rtsol];
sol[t0_?NumericQ, y0_?NumericQ, yp0_?NumericQ, w0_: 1] :=
First @ NDSolve[{ode == 0 /. w -> w0, y[t0] == y0, y'[t0] == yp0},
y, {t, 10^-20, 100}];
rtIvals[t0_, k0_] := With[{w = t0 k0},
Simplify[{a Sqrt[t0], a/(2 Sqrt[t0])} /. First @ Solve[
(-y''[t] - w^2 (y[t])/(t^2) + (y[t])^(-3) /. y -> (a Sqrt[#] &)) == 0
&& a > 0 && 1 - 4 w^2 < 0 && t > 0, a, Reals],
1 - 4 w^2 < 0 && t > 0]];
rtsol[{t0_, k0_}] := sol[t0, ##, t0 k0] & @@ rtIvals[t0, k0]
Examples
Here are several solutions with w = 1. We can see the nearby solution oscillate around the red square root graph.
Plot[
Evaluate[{y[t] /. sol[1, Sqrt[2]/3^(1/4), (1 + #) 1/(Sqrt[2] 3^(1/4))] & /@
Range[-2, 2]}],
{t, 0, 50},
PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6], Hue[0., 0.6, 1]}]
Given the nature of our analysis, a log-log plot seems like a good idea. Another approach is to subtract the square root and scale by 1/Sqrt[t], and graph as a function of s (t = E^s). The period for w = is 2 Pi/Sqrt[4 w^2 - 1]=2Pi/Sqrt[3]=3.6276. Look how well it matches. Pretty neat.
LogLogPlot[
Evaluate[{y[t] /. sol[1, Sqrt[2]/3^(1/4), (1 + #) 1/(Sqrt[2] 3^(1/4))] & /@
Range[-2, 2]}],
{t, 1/50, 50},
PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6], Hue[0., 0.6, 1]}]
Plot[Evaluate[(-(4/3)^(1/4) Sqrt[E^s] + y[E^s])/Sqrt[E^s] /.
sol[1, Sqrt[2]/3^(1/4), (1 + #) 1/(Sqrt[2] 3^(1/4))] & /@
Range[-2, 2]], {s, -Log@50, Log@50},
PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6],
Hue[0., 0.6, 1]}, GridLines -> {2. Pi/Sqrt[3] {-1, 0, 1}, None}]
Here are a few square root solutions for different values of t0, with k0 = 1.
Plot[Evaluate[{y[t] /. rtsol[{#, 1}] & /@ {1, 3, 9, 27}}],
{t, 0, 50},
PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6],
Hue[0.2, 0.6, 0.6], Hue[0., 0.6, 1]}]
NDSolveworks the way you tell it to work. That's the reality. What you should do is open up the docs and Mathematica.StackExchange and read everything you find underNDSolve\Handling singularities in DE\DE\PDE. If you have any further questions feel free to ask :) $\endgroup$Analytic solution is easy to obtain: y[t] = a Sqrt[t]? If the analytical solution is easy to find, why then not use DSolve? and get a general solution, then you can find the constant of integrations later. But using your solution even using w=1, I geteq = y[t]^-3 == y[t]/t^2 + y''[t]; sol = a Sqrt[t]; sol^(-3) - w^2 sol/t^2 - D[sol, {t, 2}]which is not zero. $\endgroup$(?)button in the top right corner of the text box. Please click it and find out how to format your posts for readability and how to place the code in a code block. $\endgroup$