function ode15s_runtime

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

  % Initial height, velocity and acceleration
  h0 = 0;
  v0 = 20;
  a0 = -9.81;
  y0 = [h0; v0];
  ode_fct = @(t, y) f (t, y, a0);

  tic;
  [t, y, te, ye, ie] = ode15s (ode_fct, [0, 5], y0, 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]
  function dydt = f (t, y, a)
    dydt = [y(2); a];
  end

  function [value, isterminal, direction] = ode_events (t, y)
    value = y(1);
    isterminal = 1;
    direction = 0;
  end

end

