[R-sig-dyn-mod] problems with times parameter in Runge-Kutta solver

[Xavier RAYNAUD]

Hello, 

I'm having problems to understand how the fixed-step Runge-Kutta solver
works in deSolve (v1.12) and more precisely how the "times" parameter
should be given. 

Below is a small example that solves Malthusian growth model using
forward euler method (I know that there are better methods for such
model but this is just to illustrate my problem and, as this model is
analytically tractable, I can compare the numerical integration to the
theoretical value): 

malthus=function(t,Y,parms) {
  with(parms,{
  dY = r*Y
  return(as.list(dY))
  })
}
parms = list(r=0.05)

y0 = 1
times = 0:100
hini = 0.001

sol = ode(y0,malthus,times=times, parms=parms,hini=hini,method="euler")

When running these lines, the population size at the end of the
integration is 148.394 which is near its theoretical value
(y0*exp(0.05*100) = 148.413). The man page for the ode function states
that the "times" argument is the "time sequence for which output is
wanted", suggesting that if I'm only interested in the value at the end
of the integration, setting times = c(0,100) should be sufficient.
However, using times = c(0,100) returns a population size of 1.648,
which is wrong. I get the correct answer if I add maxsteps = 100/hini
in the ode call. 

The diagnostics() for the ode calls with and without specifying
maxsteps only differ in  the "number of steps taken" and the "number of
function evaluation" which are respectively 1e5 (as expected from
maxsteps = 100/hini) and 10001 (which corresponds to length(times)/hini
+1 as specified in the rk() help page). The value of 1.648 I get for
time t=100 when I do not specify maxsteps is in fact the value of the
integration at t=10.001, which corresponds to the integration after
10001 steps. 

Now, my questions: 
- Am I doing things right if I want to use method="euler" or "rk4" ?
(the model I'm working on is much more complicated then malthusian
growth, and I cannot save every time steps due to memory limitations). 

- Even worse, if I use times = c(0,50,100) to obtain values for t=50
and t=100, without specifying maxsteps, returns 2.117066 for t=50 and
2.117172 for t=100  (=2.117066*(1+hini*r)) which are also wrong. The
value returned for t=50 is in fact the value corresponding to t=15.001
(after 15001 steps = length(times)/hini+1) and the value returned for
t=100 is the value corresponding to t=15.002 (15002 steps), Therefore,
once maxsteps has been reached, the solver calculates the next
requested values as if they were spaced by one integration step and
produces wrong results. I don't know if this is a bug or if I am not
using the ode function as it should be. If this is a bug, it might be
be solved by changing maxsteps default value to max(times)/hini when
using a fixed steps solver. If it is not a bug, it might help the user
to make the solver throws a warning once it reaches maxsteps to
indicate that the results could be wrong. 

Thank you very much for your help,

Regards,

Xavier.