Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|cite|improve this question
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++)
        cout << y1 << endl;
share|cite|improve this answer
Thank you, 7530! – Sean Nov 29 '11 at 13:17

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.