How to move all individuals from one class to another at specific time step Boost ODE C++

Question

I am trying to write a system of differential equations that can account for fecundity and maturation occurring at specific time steps throughout the years. Currently, this is the code that I have set up. What I want the function to do is when "mat" is equal to 1, I want all individuals from the x[0]:x[2] age classes to move to the next age class to where there are no individuals left at that age class at that time step.

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,5 > state_type;
typedef boost::array< double ,4 > params_type;


std::vector<state_type> x_result; // to store results
std::vector<double> t_result; // to store time steps

bool inSet(int value, std ::vector<int> set){
  for(int i : set){
    if(value == i){
      return true;
    }
  }
  return false;
}

void rhs( const state_type &x , state_type &dxdt , const double t, double b, NumericMatrix hazT) {
  
  double hazJ = hazT(t,0);
  double hazY = hazT(t,1);
  double hazE = hazT(t,2);
  double hazR = hazT(t,3);
  double fec = b * (inSet(t, {4, 16, 28, 40, 52}) ? 1 : 0);
  double mat = inSet(t, {3, 15, 27, 39, 51}) ? 1 : 0;
  
  dxdt[0] = (fec * x[3]) - (hazJ * x[0]) - (mat * x[0]);
  dxdt[1] = (mat * x[0]) - (hazY * x[1]) - (mat * x[1]);
  dxdt[2] = (mat * x[0]) - (hazY * x[2]) - (mat * x[2]);
  dxdt[3] = (mat * x[1]) - (hazE * x[3]);
  dxdt[4] = (mat * x[2]) - (hazR * x[4]);
  
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  //Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
  if(x_result.empty() || (std::abs(t - std::round(t)) < std::abs(t_result.back() - std::round(t)))){
    x_result.push_back(x); // save the state
    t_result.push_back(t); // save the time
  }
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
  // Clear the result vectors
  x_result.clear();
  t_result.clear();
  
  state_type x = { y[0], y[1], y[2], y[3], y[4] }; // initial conditions
  integrate_const(make_controlled( 1E-2 , 1E-2 , stepper_type () ) ,
                  [=](const state_type &x, state_type &dxdt, const double t){rhs(x, dxdt, t, b, hazT);} , 
                  x , times[0] , times[times.size() - 1] , step, write_cout);
  
  // Create a NumericMatrix to hold both time and x_values
  NumericMatrix result(x_result.size(), 6); // 3 columns for x_values and 1 column for time
  
  // Fill the NumericMatrix with time and x_values
  for (size_t i = 0; i < x_result.size(); ++i) {
    result(i, 0) = t_result[i]; // time
    result(i, 1) = x_result[i][0]; // x1
    result(i, 2) = x_result[i][1]; // x2
    result(i, 3) = x_result[i][2]; // x3
    result(i, 4) = x_result[i][3]; // x3
    result(i, 5) = x_result[i][4];
    
  }
  return result;
}

Regardless of my multiple attempts to change the function around, I cannot seem to get that behavior. The x[0] class declines but not to 0 while the x[1] and x[2] class appear to not change at all. Could someone help me understand how this function may need to be altered?

Attempted to alter the function by setting an if statement in the function; however, it did not seem to change the behavior of the function at all.

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,5 > state_type;
typedef boost::array< double ,4 > params_type;


std::vector<state_type> x_result; // to store results
std::vector<double> t_result; // to store time steps

bool inSet(int value, std ::vector<int> set){
  for(int i : set){
    if(value == i){
      return true;
    }
  }
  return false;
}

void rhs( const state_type &x , state_type &dxdt , const double t, double b, NumericMatrix hazT) {
  
  double hazJ = hazT(t,0);
  double hazY = hazT(t,1);
  double hazE = hazT(t,2);
  double hazR = hazT(t,3);
  double fec = b * (inSet(t, {4, 16, 28, 40, 52}) ? 1 : 0);
  double mat = inSet(t, {3, 15, 27, 39, 51}) ? 1 : 0;


  // Changes inserted here to move all individuals
  if(inSet(t, {3, 15, 27, 39, 51}){
  x[0] -= x[0]
  x[1] += x[0] * 0.5 - x[1]
  x[2] += x[0] * 0.5 - x[2]
  x[3] += x[1]
  x[4] += x[2]
  }
  
  // ODE updated accordingly here
  dxdt[0] = (fec * x[3]) - (hazJ * x[0]);
  dxdt[1] = (hazY * x[1]);
  dxdt[2] = (hazY * x[2]);
  dxdt[3] = (hazE * x[3]);
  dxdt[4] = (hazR * x[4]);
  
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  //Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
  if(x_result.empty() || (std::abs(t - std::round(t)) < std::abs(t_result.back() - std::round(t)))){
    x_result.push_back(x); // save the state
    t_result.push_back(t); // save the time
  }
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
  // Clear the result vectors
  x_result.clear();
  t_result.clear();
  
  state_type x = { y[0], y[1], y[2], y[3], y[4] }; // initial conditions
  integrate_const(make_controlled( 1E-2 , 1E-2 , stepper_type () ) ,
                  [=](const state_type &x, state_type &dxdt, const double t){rhs(x, dxdt, t, b, hazT);} , 
                  x , times[0] , times[times.size() - 1] , step, write_cout);
  
  // Create a NumericMatrix to hold both time and x_values
  NumericMatrix result(x_result.size(), 6); // 3 columns for x_values and 1 column for time
  
  // Fill the NumericMatrix with time and x_values
  for (size_t i = 0; i < x_result.size(); ++i) {
    result(i, 0) = t_result[i]; // time
    result(i, 1) = x_result[i][0]; // x1
    result(i, 2) = x_result[i][1]; // x2
    result(i, 3) = x_result[i][2]; // x3
    result(i, 4) = x_result[i][3]; // x3
    result(i, 5) = x_result[i][4];
    
  }
  return result;
}

Also, I am really new to c++ coding so I'm sure my methods are very rudimentary and could be streamlined in some ways.

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,5 > state_type;
typedef boost::array< double ,4 > params_type;


std::vector<state_type> x_result; // to store results
std::vector<double> t_result; // to store time steps

bool inSet(int value, std ::vector<int> set){
  for(int i : set){
    if(value == i){
      return true;
    }
  }
  return false;
}

void rhs( const state_type &x , state_type &dxdt , const double t, double b, NumericMatrix hazT) {
  
  double hazJ = hazT(t,0);
  double hazY = hazT(t,1);
  double hazE = hazT(t,2);
  double hazR = hazT(t,3);
  double fec = b * (inSet(t, {4, 16, 28, 40, 52}) ? 1 : 0);
  double mat = inSet(t, {3, 15, 27, 39, 51}) ? 1 : 0;
  
  dxdt[0] = (fec * x[3]) - (hazJ * x[0]) - (mat * x[0]);
  dxdt[1] = (mat * x[0]) - (hazY * x[1]) - (mat * x[1]);
  dxdt[2] = (mat * x[0]) - (hazY * x[2]) - (mat * x[2]);
  dxdt[3] = (mat * x[1]) - (hazE * x[3]);
  dxdt[4] = (mat * x[2]) - (hazR * x[4]);
  
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  //Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
  if(x_result.empty() || (std::abs(t - std::round(t)) < std::abs(t_result.back() - std::round(t)))){
    x_result.push_back(x); // save the state
    t_result.push_back(t); // save the time
  }
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
NumericMatrix boostExample_t(NumericVector y, NumericVector times, double step, double b, NumericMatrix hazT) {
  // Clear the result vectors
  x_result.clear();
  t_result.clear();
  
  state_type x = { y[0], y[1], y[2], y[3], y[4] }; // initial conditions
  integrate_const(make_controlled( 1E-2 , 1E-2 , stepper_type () ) ,
                  [=](const state_type &x, state_type &dxdt, const double t){rhs(x, dxdt, t, b, hazT);} , 
                  x , times[0] , times[times.size() - 1] , step, write_cout);
  
  // Create a NumericMatrix to hold both time and x_values
  NumericMatrix result(x_result.size(), 6); // 3 columns for x_values and 1 column for time
  
  // Fill the NumericMatrix with time and x_values
  for (size_t i = 0; i < x_result.size(); ++i) {
    result(i, 0) = t_result[i]; // time
    result(i, 1) = x_result[i][0]; // x1
    result(i, 2) = x_result[i][1]; // x2
    result(i, 3) = x_result[i][2]; // x3
    result(i, 4) = x_result[i][3]; // x3
    result(i, 5) = x_result[i][4];
    
  }
  return result;
}