Listing 5: Core functions of the Gear integrator.

template <class C, class T, class B> void
Gear<C,T,B>::step(void)
{
    // Iterate until xnext is self-consistent.
    xnext = statevector;
    unsigned int iter = 0;
    do {
        xtemp = xnext;
        return_derivs(xnext,t+h,dxdt);
        (this->*(calc_xnext[curr_order-1]))();
        iter++;
        if (iter > maxiter)
          throw std::runtime_error
            ("Too many iterations in Gear method.");
    } while (!approx_equal(xtemp,xnext));

    // Past state queue management
    past_state.push_back(statevector);
    if (order > curr_order) curr_order++;
    else past_state.pop_front();

    // Update current state
    statevector = xnext;
    t += h;
}

// First order Gear formula: xnext = statevector + h*dxdt
template <class C, class T, class B>
void Gear<C,T,B>::Gear1xnext(void)
{
    xnext = dxdt;    xnext *= B(h);    xnext += statevector;
}

// Second order Gear formula:
// xnext = (4*statevector - past_state[0] + 2*h*dxdt)/3
template <class C, class T, class B>
void Gear<C,T,B>::Gear2xnext(void)
{
    xnext = dxdt;    xnext *= B(2.0*h);
    xnext += B(4.0)*statevector;    xnext -= past_state[0];    
    xnext /= B(3.0);
}

// Gear3xnext and Gear4xnext not shown

template <class C, class T, class B>
void Gear<C,T,B>::set_step_size(T new_h)
{
    past_state.clear();
    curr_order = 1;
    h = new_h;
}