/* Output: f is the i-th difference of the curve after one step  */
static void fd_fwd(double f[4])
/* Output: f is the i-th difference of the curve after one step  */
static void fd_fwd(double f[4])
{
    f[0] += f[1];
    f[1] += f[2];
    f[2] += f[3];
}