Hello, Question in brief: what's a good way, using TS, to solve a system of ODEs which looks like this? M U_t = -K U + f(t) In addition, how can the code be kept nearly identical to allow different TS types to be used (EULER, BEULER, TSRUNGE_KUTTA, TSCRANK_NICHOLSON, etc)? I'd also like to take possible advantages of M and K being constants. Background and what I've tried: If f(t) were not present TSSetMatrices seems like the convenient way ahead. In fact, I can still use the same function in presence of f(t) for an explicit time-stepping method like RK or Euler if I use a MatShell for Arhs, and sneak in f(t) there. But I cannot use an implicit method like Backward Euler with a MatShell for Arhs easily. It needs MatFactor to be implemented for the shell matrix (at least). Is there something I'm missing? This is a common system of equation. Most TS examples solve u_t = u_xx, convenient for TSSetMatrices. ts/examples/tests/ex3.c looks like it should work for a system with f(t). But it needs the Jacobian matrix of RHS and it evaluates it using finite differences. In addition, for each time-step it computes an ILU factorization. Thanks, Chetan