# 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...

-
So, MATLAB. It'd be nice if you'd mentioned the MATLAB code you already have. – J. M. Nov 29 '11 at 5:05
Actually no. I used R's solver. Basically the function is defined as $$dg = y[2]; dg^{2} = \frac{a(y[2])^2}{by[2]+cy[1]+d}.$$ – Sean Nov 29 '11 at 5:16
Fine, R. But no code? – J. M. Nov 29 '11 at 5:17

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;
}
}

-
Thank you, 7530! – Sean Nov 29 '11 at 13:17