Boost OdeInt cannot set class property when entering integration routine

I'm using boost::numeric::odeint to build a custom integrator. Other than computing the discrete derivative of the state, I'm computing the output equation during each integration step entering operator (). I do know it's not efficient, but as of now I need it to understand how the library is working.

While debugging, inside the operator () method both y and test variables are computed fine. But, as soon as the program returns to main loop, sys objects has y and test variables uninitialized.

Below the code is presented and the result that reproduces what discussed.

TestSystem.h

#pragma once

#include <Eigen/Dense>
#include <iostream>

using namespace std;
using namespace Eigen;

const int n = 2; // State size
const int q = 1; // Output size
const int m = 1; // Input size

typedef Vector<double, n> state_type;
typedef Vector<double, m> input_type;
typedef Vector<double, q> output_type;

typedef Matrix<double, n, n> state_matrix_type;
typedef Matrix<double, n, m> input_matrix_type;
typedef Matrix<double, q, n> output_matrix_type;
typedef Matrix<double, q, m> feedthrough_matrix_type;

class TestSystem
{

public:
    state_matrix_type A{ 
        { -4.0, -3.0}, 
        { 1.0, 0.0} 
    };
    input_matrix_type B{ 1.0, 0.0 };
    output_matrix_type C{ 1.0, 1.0 };
    feedthrough_matrix_type D{ 0.0 };

    double test;

    input_type u;
    output_type y;
    void operator() (const state_type& x, state_type& x_dot, const double t);
};

TestSystem.cpp

#include "TestSystem.h"

void TestSystem::operator()(const state_type& x, state_type& x_dot, const double t)
{
    x_dot = A * x + B * u;

    this->y = C * x + D * u;

    this->test = (C * x + D * u)(0);

    cout << "Output computed within integration step: " << endl;

}

main.cpp

#include <iostream>
#include <Eigen/Dense>
#include "TestSystem.h"

using namespace Eigen;
using namespace std;
using namespace cnpy;

const double PI = 3.14159265359;

const int N = 5000;

VectorXd u(N);

#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>

using namespace boost::numeric::odeint;

int main()
{
    VectorXd t(N);

    t = VectorXd::LinSpaced(N, 0.0, 1.0);

    TestSystem sys;

    u = (2 * PI * 200 * t).array().cos() + (2 * PI * 20 * t).array().sin();

    state_type x = state_type::Zero();
    
    // Just one step for example sake
    for (size_t i = 0; i < 1; i++) {
        input_type _u{ u(i) };
        // _u << u(i);

        sys.u = _u;

        runge_kutta4<state_type> rk;

        rk.do_step(sys, x, t(i), t(1));

        cout << "System object output variable: " << sys.test << endl;
    }
}

Result

Output computed within integration step: 0
Output computed within integration step: 0.00010002
Output computed within integration step: 9.999e-05
Output computed within integration step: 0.00019998
System object output variable: -9.25596e+61

Solution

The system is passed by value. To force reference semantics, use a reference wrapper:

Live On Compiler Explorer

#include <Eigen/Dense>
#include <iostream>

using Eigen::Matrix;
using Eigen::Vector;
using Eigen::VectorXd;

int constexpr n = 2; // State size
int constexpr q = 1; // Output size
int constexpr m = 1; // Input size

using state_type  = Vector<double, n>;
using input_type  = Vector<double, m>;
using output_type = Vector<double, q>;

using state_matrix_type       = Matrix<double, n, n>;
using input_matrix_type       = Matrix<double, n, m>;
using output_matrix_type      = Matrix<double, q, n>;
using feedthrough_matrix_type = Matrix<double, q, m>;

struct TestSystem {
    state_matrix_type       A{{-4.0, -3.0}, {1.0, 0.0}};
    input_matrix_type       B{1.0, 0.0};
    output_matrix_type      C{1.0, 1.0};
    feedthrough_matrix_type D{0.0};

    double test;

    input_type  u;
    output_type y;
    void        operator()(state_type const& x, state_type& x_dot, double t);
};

void TestSystem::operator()(state_type const& x, state_type& x_dot, double) {
    x_dot = A * x + B * u;

    this->y = C * x + D * u;

    this->test = (C * x + D * u)(0);
    std::cout << "Output computed within integration step: " << test << std::endl;
}

#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>

namespace ode = boost::numeric::odeint;

int main() {
    int constexpr N = 5000;
    VectorXd t(N);

    t = VectorXd::LinSpaced(N, 0.0, 1.0);

    TestSystem sys;

    VectorXd u(N);
    u = (2 * M_PI * 200 * t).array().cos() + (2 * M_PI * 20 * t).array().sin();

    state_type x = state_type::Zero();
    
    // Just one step for example sake
    for (size_t i = 0; i < 1; i++) {
        input_type _u{ u(i) };
        // _u << u(i);

        sys.u = _u;

        ode::runge_kutta4<state_type> rk;

        rk.do_step(std::ref(sys), x, t(i), t(1));

        std::cout << "System object output variable: " << sys.test << std::endl;
    }
}

Printing

Output computed within integration step: 0
Output computed within integration step: 0.00010002
Output computed within integration step: 9.999e-05
Output computed within integration step: 0.00019998
System object output variable: 0.00019998