Runge-Kutta fixed step solvers

This function implements a fixed-step Runge-Kutta solver for explicit and implicit methods (and with optional adaptive step size control).

The function supports both explicit and implicit methods, and embedded methods as well. Any Runge-Kutta method can be added simply by specifying their butcher tableaus. The algorithm itself is generic and relatively compact. About 34 methods are currently implemented.

MATLAB's ODE solvers are all variable-step and don't even offer an option to run with fixed step size. This is because adaptive step size can make a solver both faster and more precise compared to fixed step size. However, sometimes there are good reasons to choose a fixed-step solver:
- parameter studies (comparing simulation results for different model parameters)
- calculating finite-difference jacobians of the simulation results (adaptive step size control can introduce significant noise)
- performing point-wise calculations where solver output and measurement data must refer to the same time vector
- having pre-allocated arrays for simulation results and fixed computing time

Interface and options are explained in comments. There are two examples:
example1 solves a damped and driven harmonic oscillator with different methods and step sizes for comparison of performance.
example2 solves the duffing equation with two different implicit solvers.

Example (harmonic oscillator):
t=linspace(0,10);
y=rkfs(@(~,x)[x(2) -x(1)],t,[1 0]);
plot(t',y');