In the steppers, there are currently routines to compute the residual and the Jacobian. We should add an "applyJacobian" routine to enable, e.g., matrix-free solvers. This could be very important if we want to use Pressio to time march certain types of FOMs. This could look like,
void applyJacobian(vector_t v , vector_t Jv){
/* compute Jv through appropriate method, e.g., for finite difference
Jv = 1./eps*( residual(v + eps) - residual(v) )
*/
}