In what way depend the results on the step size? In you example, the step size that you pass to the integrate function is in principle the observation time step. But it is also the initial time step which might be in appropriate for your use case.
After running a few more tests. It seems like the solution does not vary when step_size is changed but sometimes the solver will freeze, as I described earlier. For example, with my original code that uses the runge_kutta_fehlberg78 the solver gives the same answer when stepsize < 0.7 but when it it larger it will freeze and give no answer at all.
Can you try to solve the above example with integrate times:
vector< double > obs_times = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 }; double step_size = 0.02; // now it is really the initial step size integrate_times( make_controlled< stepper >( abs_err , rel_err ) , model2 , x , step_size , times.begin() , times.end() , Print_observer );
I used integrate_times as you described above but it also froze after
the .1 step. After having waited for a very long time I eventually get
this error message
"""
terminate called after throwing an instance of 'boost::exception_detail::clone_impl
If you system is stiff you should also look at the Rosenbrock method. It is much better suited for stiff problems. But its needs to calculate the Jacobian.
The biggest issue is that the Jacobian is too complicated to calculate, it depends on some population dynamics and higher order interpolation. Is there no other way to solve stiff systems. I've read about other methods such as CVODE and LSODA, but couldn't figure out how to get those to work.