Numerically solve second-order ODE I want to solve a second-order ODE in the form of
$$ y^{''} = \frac{a (y^{'})^2}{b y^{'}+cy+d} $$ by numerical method (eg, solver ODE45), given initial condition of $y(0)$ and $y'(0)$. The results are wield and numbers go out of machinery bound. I guess the catch is that what is in the denominator becomes highly unstable when it converges to zero. I tried bound it away from zero with no avail.
Could anyone provide insights on how to proceed with the numerical procedure? Thanks in advance...

Discretizing the ODE by finite differences gives
$$\frac{y_2-2y_1 + y_0}{h^2} = \frac{a\left(\frac{y_1-y_0}{h}\right)^2}{b\left(\frac{y_1-y_0}{h}\right) + cy_1 + d},$$
or
$$y_2 = 2y_1 - y_0 + \frac{ah(y_1-y_0)^2}{b(y_1-y_0)+chy_1+dh}.$$
Here's C++ code I wrote which has no trouble integrating this ODE for $a=-10,b=c=d=1$, initial conditions $y(0)=0,y'(0)=10$ and time step $h=0.01$. I'm sure you can adapt it to whatever language you prefer:
#include <iostream>

using namespace std;

double y0, y1;

void step(double dt)
{
    double y2 = 2*y1-y0 - 10*dt*(y1-y0)*(y1-y0)/(y1-y0+dt*y1+dt);
    y0 = y1;
    y1 = y2;
}

int main()
{
    y0 = 0;
    y1 = .1;
    for(int i=0; i<1000; i++)
    {
        step(0.01);
        cout << y1 << endl;
    }
}

