From owner-bug-octave at bevo dot che dot wisc dot edu  Mon Oct 28 15:10:33 1996
Subject: Bug in Octave's LSODE-function
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Janne Karonen <janne dot karonen at kemira dot com>
cc: bug-octave at bevo dot che dot wisc dot edu
Date: Mon, 28 Oct 1996 15:10:20 -0600

On 17-Oct-1996, Janne Karonen <janne dot karonen at kemira dot com> wrote:

: To recreate the bogus behaviour issue the following commands in Octave's
: command window. Before that make sure you have m-file bug.m in the
: working directory.
: 
: 	lsode_options("relative tolerance",1e-3)
: 	lsode_options("absolute tolerance",1e-3)
: 
: 	data=lsode("bug",[0],(t=linspace(0,20,21)'));
: 
: This will lead to the following error messages.
: 
: 	error: matrix index = 23 exceeds maximum dimension = 21
: 	error: evaluating index expression near line 9, column 3
: 	error: evaluating assignment expression near line 9, column 2
: 	error: called from 'bug' in file '/home/jkaronen/octave/bug.m'
: 	error: lsode: evaluation of user-supplied function failed
: 
: Simulation of the 'bug.m' can be performed if the maximum step size
: is set to 1 with the following command.
: 
: 	lsode_options("maximum step size",1)
: 
: When 'bug.m' file is modified commenting out the line which returns
: the simulation time 't' on every round it can be seen that the last
: value of 't' is 22,289. This is for some unknown reason over the end
: point of the simulation.

This is a feature of the way LSODE is implemented.  It may choose
points beyond the last value of tout specified and interpolate for the
values to return.

The Fortran subroutine LSODE does allow you to tell it to not
integrate past a specified point but Octave doesn't currently have a
way for you to specify the stopping point, though I've added this to
the PROJECTS file.

: ---------------------------------------------------------------------
: function f=bug(x,t)
: 
: a=[5 5 5 6 6 6 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]';
: 
: # Comment this out to see how Octave chooses simulation time steps.
: # t
: 
: index=(t-(t-floor(t)))+1;
: 
: b=a(index);
: 
: f=1/(1+x^2)-b*x^2;
: 
: endfunction
: ---------------------------------------------------------------------


I would recommend modifying your function to calculate the index like
this

  index=min((t-(t-floor(t)))+1, length(a));

instead.

jwe