RFC: solver ode23s for stiff problems

The implementation is based on the "modified Rosenbrock triple" described in
> L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22

which is also the basis for MATLAB's `ode23s`. In addition to the usual keywords, the solver accepts a keyword `jacobian`, which can be used to supply a function giving the Jacobian of the RHS. Otherwise a (very crude) finite-difference method is used (pulled out from `oderosenbrock`).

A simple example is:
```
# Example: van der Pol oscillator
import ODE: ode23s

function vdP(t,y,mu)
    dy = zeros(y)

    dy[1] = y[2]
    dy[2] = -y[1] + mu*(1 - y[1]^2)*y[2]

    return dy
end

function jac_vdP(t,y,mu)
    dFdy = zeros(length(y), length(y))

    dFdy[1,1] = 0.
    dFdy[1,2] = 1.

    dFdy[2,1] = -1. - mu*2*y[1]*y[2]
    dFdy[2,2] = mu*(1 - y[1]^2)

    return dFdy
end

tout, yout = ode23s((t,y)->vdP(t,y,1000), [2.,0.], [0., 3000.], jacobian=(t,y)->jac_vdP(t,y,1000), reltol=1e-3, abstol=1e-6)
```
Yielding for the first component of `yout`
![ode23s-vdp1000](https://cloud.githubusercontent.com/assets/6170781/4168694/96fff46c-351f-11e4-98bb-9545faa66977.png)

This PR should also "fix" #1.