Feature: Add possibility to call direct implementation of RHS

[ranocha]

It would be nice to add the option to provide a direct implementation of the standard ODE RHS for production-destruction systems. This will typically be more efficient than the automatic conversion from the PDS system.

This will make it easier to implement fair comparisons with standard methods just calling the ODE RHS and also improve the performance of methods that use the current version of dense output (that should probably be changed anyway, see #52).

[ranocha]

julia> struct Foo{P}
           p::P
       end

julia> (foo::Foo)(u) = sum(foo.p(u))

julia> foo = Foo(u -> ones(2))
Foo{var"#1#2"}(var"#1#2"())

julia> foo([1, 2])
2.0

julia> 

julia> 

julia> struct Bar{P, RHS}
           p::P
           rhs::RHS
       end

julia> Bar(p) = Bar(p, rhs_default)
Bar

julia> (bar::Bar)(u) = bar.rhs(u)

julia> Bar(p) = Bar(p, u -> rhs_default(p, u))
Bar

julia> rhs_default(p, u) = sum(p(u))
rhs_default (generic function with 1 method)

julia> bar = Bar(u -> ones(2))
Bar{var"#7#8", var"#5#6"{var"#7#8"}}(var"#7#8"(), var"#5#6"{var"#7#8"}(var"#7#8"()))

julia> bar([1, 2])
2.0

julia> bar2 = Bar(u -> ones(2), u -> 2)
Bar{var"#9#11", var"#10#12"}(var"#9#11"(), var"#10#12"())

julia> bar([1, 2])
2.0

julia> bar2([1, 2])
2