x0 = 2; y0 = 1;  % initial values
t = [0,linspace(6,7,1000)];  % examine times between 6 and 7

lsode_options("relative tolerance",1e-10)
lsode_options("absolute tolerance",1e-10)
XY = lsode('VolterraLotkaC',[x0,y0],t); % solve the ODE

y = XY(2:end,1)-x0; t = t(2:end);      % examine the first component only
plot(t,y);  % visual test for zero
grid on

s = sign((XY(2:end,1)-x0).*(XY(1:end-1,1)-x0)); % detect sign changes
pos = find(s<0);                              % position of sign change
% use linear interpolation to determine the partial time step
dt = (t(pos)-t(pos-1))*y(pos-1)/(y(pos-1)-y(pos))
T  = t(pos-1)+dt                             % estimate the period
XYt = lsode('VolterraLotkaC',XY(pos,:),[T-dt,T]);  % compute the value at T
XYt(2,:)-[x0,y0]    % difference with initial values


