Exa = 1;

switch Exa
case 1   % use ode45 from Matlab or Octave
%  function dy = ODEPend(t,y)
%    k  = 1; alpha = 0.1;
%    dy = [y(2);-k*y(1)-alpha*y(2)];
%  endfunction
  vopt  = odeset('RelTol', 1e-10, 'AbsTol', 1e-10,'NormControl', 'on');
  [t,y] = ode45(@ODEPend, [0 25], [0 0.9], vopt);
  figure(1); plot(t,y);
  xlabel('time');ylabel('position and velocity');grid on
  figure(2); plot(y(:,1),y(:,2));
  xlabel('position'); ylabel('velocity');  grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 2   % use the function rk45.m, working with Octave or Matlab
  [t,y] = rk45('ODEPend',0,25,[0 0.9],1e-5,1e-5);
  figure(1); plot(t,y);
  xlabel('time');ylabel('position and velocity');grid on
  figure(2); plot(y(:,1),y(:,2));
  xlabel('position'); ylabel('velocity');  grid on
end
