So, it turns out that scipy doesn't really implement its own ode solver, but instead wraps a bunch of FORTRAN code, and just silently drops like 80% of the functionality of that code without documenting much of what it kept and what it ignored. Anyway, after sifting through a bunch of forums and a whole bunch of source code, I came to the conclusion that the only function that can actually do adaptive time stepping and supply dense output is the wrapper around the vode integrator. The other thing that I needed was to detect when my trajectory crosses a certain boundary defined by $\phi(y) = 0$ and trigger an event at that point. In my case I simply needed to restart integration after adding the crossing point to the solution array. The following code does what I needed:
def f(t, y):
# define your \dot{y}=f(y,t) here
def phi(y):
# this function defines the boundary that is crossed
from scipy.integrate import ode
# set up the solver, with a maximum time step of 1e-1
# an initial condition of y0 = init and t0 = 0
solver = ode(f)
solver.set_integrator('vode', max_step=1e-1)
solver.set_initial_value(init, 0)
tfinal = 30 # run this for 30 seconds
t = [0]
y = [init]
# call the solver iteratively with the 'step' option
# what this does is it returns the result of the integration
# after the next (variable) time step, and the length of that step
while solver.successful() and solver.t < tfinal :
solver.integrate(tfinal, step=True)
y.append(solver.y);
t.append(solver.t);
# check for zero crossing
gg1 = phi(y[-1])
gg0 = phi(y[-2])
if gg0*gg1 < 0 :
# find zero crossing
# simple linear interpolation will do for now
tcross = t[-2] - gg0*(t[-1]-t[-2])/(gg1-gg0)
ycross = y[-2] - gg0*(y[-1]-y[-2])/(gg1-gg0)
# replace wrong y and t
y[-1] = ycross
t[-1] = tcross
# reset integration and start again
solver.set_initial_value(ycross, tcross)
print("crossed the boundary at time ", tcross)
