Imaging, you want to numerically integrate a harmonic oscillator with friction. The equations of motion are given by x'' = -x + γ x'. Odeint only deals with first order ODEs that have no higher derivatives than x' involved. However, any higher order ODE can be transformed to a system of first order ODEs by introducing the new variables q=x and p=x' such that w=(q,p). To apply numerical integration one first has to design the right hand side of the equation w' = f(w) = (p,-q+γ p):

/* The type of container used to hold the state vector */
typedef std::vector< double > state_type;

const double gam = 0.15;

/* The rhs of x' = f(x) */
void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
{
    dxdt[0] = x[1];
    dxdt[1] = -x[0] - gam*x[1];
}

Here we chose vector<double> as the state type, but others are also possible, for example boost::array<double,2>. odeint is designed in such a way that you can easily use your own state types. Next, the ODE is defined which is in this case a simple function calculating f(x). The parameter signature of this function is crucial: the integration methods will always call them in the form f(x, dxdt, t) (there are exceptions for some special routines). So, even if there is no explicit time dependence, one has to define t as a function parameter.

Now, we have to define the initial state from which the integration should start:

state_type x(2);
x[0] = 1.0; // start at x=1.0, p=0.0
x[1] = 0.0;

For the integration itself we'll use the integrate function, which is a convenient way to get quick results. It is based on the error-controlled runge_kutta54_cash_karp stepper (5th order) and uses adaptive step-size.

size_t steps = integrate( harmonic_oscillator ,
        x , 0.0 , 10.0 , 0.1 );

The integrate function expects as parameters the rhs of the ode as defined above, the initial state x, the start-and end-time of the integration as well as the initial time step=size. Note, that integrate uses an adaptive step-size during the integration steps so the time points will not be equally spaced. The integration returns the number of steps that were applied and updates x which is set to the approximate solution of the ODE at the end of integration.

It is also possible to represent the ode system as a class. The rhs must then be implemented as a functor - a class with an overloaded function call operator:

/* The rhs of x' = f(x) defined as a class */
class harm_osc {

    double m_gam;

public:
    harm_osc( double gam ) : m_gam(gam) { }

    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
    {
        dxdt[0] = x[1];
        dxdt[1] = -x[0] - m_gam*x[1];
    }
};

which can be used via

harm_osc ho(0.15);
steps = integrate( ho ,
        x , 0.0 , 10.0 , 0.1 );

In order to observe the solution during the integration steps all you have to do is to provide a reasonable observer. An example is

struct push_back_state_and_time
{
    std::vector< state_type >& m_states;
    std::vector< double >& m_times;

    push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
    : m_states( states ) , m_times( times ) { }

    void operator()( const state_type &x , double t )
    {
        m_states.push_back( x );
        m_times.push_back( t );
    }
};

which stores the intermediate steps in a container. Note, the argument structure of the ()-operator: odeint calls the observer exactly in this way, providing the current state and time. Now, you only have to pass this container to the integration function:

vector<state_type> x_vec;
vector<double> times;

steps = integrate( harmonic_oscillator ,
        x , 0.0 , 10.0 , 0.1 ,
        push_back_state_and_time( x_vec , times ) );

/* output */
for( size_t i=0; i<=steps; i++ )
{
    cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
}

That is all. You can use functional libraries like Boost.Lambda or Boost.Phoenix to ease the creation of observer functions.

The full cpp file for this example can be found here: harmonic_oscillator.cpp