function ode15sfun


  opt = odeset("events", @ode_events, "abstol", 1e-6, "reltol", 1e-6, "maxstep", 0.01);
  % [t,y,te,ye,ie] = ode15s (model, tspan, y0, options)
  [t, y, te, ye, ie] = ode15s (@f, [0, 5], [0; 20], opt);

  figure;
  plot (t, y(:,1), ".");
  xlabel ("t");
  ylabel ("y");
  grid on;

  % y'' = a
  % Let y1 := y and y2 := y' then
  % [y1', y2'] = [y2, a]
  function dydt = f(t, y)
    dydt = [y(2); -9.8];
  end

  % ode15s event function has to input arguments
  function [value, isterminal, direction] = ode_events (t, y)
    value = y(1);
    isterminal = 1;
    direction = 0;
  end

end

