Describe the bug

The integral collocation method implementation of this paper is numerically unstable.
We use the Chebyshev transform class provided by Boost to compute the fixed point operator (Lagrange interpolation at Chebyshev nodes).

The implementation of the ODE solver is located at include/ode_solvers/integral_collocation.hpp.

To Reproduce
Steps to reproduce the behavior:

Simulate the ode x''(t) = -x(t) with x(0) = 0 and x'(0) = 1.
The numerical solver returns a sine function with increasing amplitude.
I believe this has to be with the parameterization of the algorithm.

Expected behavior

The solution is x(t) = sin(t) which is bounded above by 1.