function ode15ifun

  opt = odeset ("events", @ode_events, "abstol", 1e-6, "reltol", 1e-6, ...
                "maxstep", 0.00004);

  % Initial height, velocity and acceleration
  h0 = 0;
  v0 = 20;
  a0 = -9.81;
  t0 = 0;
  y0 = [h0; v0];
  yp0 = [v0; a0];
  % Compute consistent initial values
  ode_fct = @(t, y, ydot) f (t, y, ydot, a0);
  [y0_new, yp0_new] = decic (ode_fct, t0, y0, [1; 1], yp0, [0; 0]);

  tic;
  [t, y, te, ye, ie] = ode15i (ode_fct, [0, 5], y0_new, yp0_new, opt);
  toc;

  length (y)
  figure;
  hold on;
  yexct = 0.5 * a0 * t.^2 + v0 * t + h0;
  plot (t, yexct, "r.");
  plot (t, y(:,1), "b.");
  xlabel ("t");
  ylabel ("y");
  legend ({"exact", "ode15"});
  title ("Ball");
  grid on;

  % ODE: y'' = a
  % Let y1 := y and y2 := y' then
  % [y1', y2'] = [y2, a]
  % Rearranging gives
  % 0 = [y2 - y1', a - y2']
  function dydt = f (t, y, ydot, a)
    dydt = [y(2) - ydot(1); a - ydot(2)];
  end

  function [value, isterminal, direction] = ode_events (t, y, yp)
    %fprintf ("t: %f\th: %f\tv(y): %f\tv(yp): %f\ta: %1.2f\n", t, y(1), y(2), yp(1), yp(2));
    value = y(1);
    isterminal = 1;
    direction = 0;
  end

end

