function ode_2a  % ode_2a.m
  options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5]);
  [T,Y] = ode45(@rhs, [0 4], [0 0 0], options);
  plot(T, Y(:,1), '-', T, Y(:,2), '-.', T, Y(:,3), '.')
  title('ode_2a.m');

  function dy=rhs(t, y)
   dy = zeros(3,1); % force a column vector
   dy(1) = 0.5;  % dv/dt = acceleration, y(1) is velocity
   dy(2) = y(1); % dx/dt = v,            y(2) is position
   dy(3) = y(2); % ds/dt = x,            y(3) is integral of position
  end % rhs
end % ode_2a
