I want to know what is the correct way of solving an ODE of the following form using `DifferentialEquations`

.

Consider for example the ODE satisfied by a function u(r) (prime indicates a derivative)

Such ODEs often arise when solving problems under the assumption of spherical or cylindrical symmetry. Solutions are found by assuming a power series form for u. We can verify that the function u(r) = gr satisfies the ODE as well as the initial conditions.

A typical way of solving such a second order ODE numerically is to re-write it as a system of two first order ODEs by introducing a variable v = u'. The system is then:

Using `DifferentialEquations`

I attempt to solve the above problem as follows:

```
using DifferentialEquations
function g(du, u, p, t)
du[1] = u[2]
du[2] = -u[2]/t + u[1]/t^2
end
u0 = [0.0, -1.0]
tspan = (0.0, 1.0)
prob = ODEProblem(g, u0, tspan)
sol = solve(prob)
```

This produces the following warnings,

```
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uJP6R/src/initdt.jl:74
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uJP6R/src/solve.jl:376
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/EFqMn/src/integrator_interface.jl:162
```

I understand that this is probably because the differential equation cannot be evaluated at r = 0 as this will produce `NaN`

s.

I can offset `tspan`

to start a little bit away from r = 0 and produce a result. My question, however, is how does one go about solving such an ODE using `DifferentialEquations`

?

**Note**: For the function u to be analytic at r = 0 we can show by Taylor expansion that u''(0) = 0. I use this condition in my regular Finite Difference solver. Perhaps I can introduce such a constraint somehow into the `ODEProblem`

?