boost::odeint called within member class

this is a personal project I've been working on and I can't figure out what's going on here (just learning C++). I found answers to very similar problems, but I can't seem to execute the solution. Here is my code with some of the unimportant bits trimmed out:

#include <iostream>
#include <cmath>
#include <complex>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include <gsl/gsl_roots.h>

class Riemann
{
public:
    // constructor
    Riemann(double leftP, double rightP, double leftrho, double rightrho, \
        double leftvx, double rightvx, double leftvy, double rightvy, double gam);

    double PL,PR,rhoL,rhoR,vxL,vxR,vyL,vyR,gamma;

    // function prototypes
    double shockvelocity(double Pg, int sign);
    double rarefactionvelocity(double Pg, int sign);
    void RfODE(const boost::array<double,6> &vrhovt, \
        boost::array<double,6> &dvrhovtdp, double t);

//  ~Riemann();
};

Riemann::Riemann(double leftP, double rightP, double leftrho, double rightrho, \
        double leftvx, double rightvx, double leftvy, double rightvy, double gam){
    // constructs Riemann public variables
}

double Riemann::shockvelocity(double Pg,int sign){
    // calculate a shock velocity, not important here...
}



void Riemann::RfODE(const boost::array<double,6> &vrhovt, \
    boost::array<double,6> &dvrhovtdp, double t){
    // calculates the ODE I want to solve

}

double Riemann::rarefactionvelocity(double Pg, int sign){
    double dpsize=0.00001;
    double P,rho,vx,vy,vtest;
    //
    boost::array<double,6> vrhovt = {vx,rho,vy,double(sign),P,gamma}; // initial conditions
    boost::numeric::odeint::integrate(std::bind(&Riemann::RfODE,std::ref(*this),std::placeholders::_1, 
        std::placeholders::_2, std::placeholders::_3),vrhovt,P,Pg,dpsize);
    std::cout<<"vRarefaction="<<vrhovt[0]<<std::endl;
    return vrhovt[0];
}

double FRiemann(double Pg, void* Riemannvalues){
    Riemann* Rvals = (Riemann*)Riemannvalues;
    // calls on Riemann::rarefactionvelocity at some point
}


int main(){
    double PL= 1000.0;
    double PR= 0.01;
    double rhoL= 1.0;
    double rhoR= 1.0;
    double vxL= 0.0;
    double vxR= 0.0;
    double vyL= 0.0;
    double vyR= 0.0;
    double gam = 5.0/3.0;

    // calls FRiemann to get a root

}

What's happening is the code is going through, calling Riemann::rarefactionvelocity just fine, but for some reason RfODE is never executed (ex. print statements in this function never execute) and the value for vrhovt[0] returned is of course the value it began with, vx. No compiler errors, either (using gcc 4.8.1 and -std=c++11 and -O2 tags) This is very strange because I've tested the rarefaction-specific functions on their own (outside of the Riemann class) and they work -- the problem seems to be that they're in this class. Given how Riemann solvers work, though, I had my reasons for making a class out of these functions and really would like to find a way to make this work without doing a massive rewrite and changing the class structure.

Any help is much appreciated! Thank you! : )

Solution

It might be possible that P is not initialized correctly. At least I don't see it in your code. P needs to be smaller than PG otherwise your are already behind your the end of the integration.

Also, don't use bind, use a lambda instead. I think bind is somehow obsolete in C++11/C++14. It might be possible that bind doesn't get the references correct.

double Riemann::rarefactionvelocity(double Pg, int sign)
{
    // ...

    // not tested
    using namspace boost::numeric::odeint;
    integrate( [this](auto const& x, auto &dxdt ,auto t ) {
        this->RfODE(x, dt, t); } ,vrhovt,P,Pg,dpsize);
}