% This is a front end for ODE-solving using any solver that takes the same
% arguments as euler.m.  Step size is fixed.
function [x,t] = ode_fixed(solver, ode, x0, t0, t1, Nsteps)

h = (t1 - t0)/Nsteps;
t = (t0 + [0:Nsteps]*h)'; % column vector of times
Ndim = numel(x0); % How big is the phase space?
x = zeros(Nsteps+1, Ndim); % each row is a state at a given timestep.

x(1,:) = x0;

for n = 2:Nsteps+1
   x(n,:) = solver(ode, x(n-1,:), t(n-1), h);
      % This is counterintuitive --- 
      % pass function handle as ode rather than %ode.
end