function [ y,t ] = ode(D,y_0,t_0,t_f,N_s,ord,fig) 

% time stepping, to solve the D-vector ODE dy/dt = f(y,t),
%    the function f(y,t) must be provided as f_RHS.m

% INPUT
% dimension: D
% initial data: y_0
% t_0:  initial  time
% t_f:  final time  
% N_s: number of time steps
% ord: order of integrator (1=Euler,2=Heun,4=Runge-Kutta)
% fig: set equal to 1 for figures

% OUTPUT
% grid of times: t (N_s x 1 vector)
% solution history: y (N_s x D matrix)

% the function returns the values of y and t in a file 'data' 
% the function also plots the components of the solution y  vs. t 
  % the function final prints t,y to a file 'data'

dt = (t_f - t_0)/N_s; % size of time step
t = t_0:dt:t_f; % the values of t for which solution is required 

% initialise y to be same size as t 
N = length(t); % note! N=N_s+1
y = zeros(N,D);

j = 1; % j will be the time index (1 <= j <= N)
y(j,:) = y_0(:)'; % initial condition for first value of t 

while t(j) < t_f 
% code to advance solution one step 
% using estimate of average value of dxdt over interval 
% details of average value to be worked out in function dydt
yj=y(j,:);
y(j+1,:) = yj  + dt*dydt(ord,yj,t(j),dt); 
      % NOTE: average rate dydt may depend on t, y and dt
j = j + 1;
end

if fig==1

for k=1:D,
figure
k
z=y(:,k);
plot(t,z) % the 1st component of solution 
xlabel('time t')
%ylabel('y_k') 
ylabel(sprintf('y_%d', k))
pause
end

if D==2
figure
k
z=y(:,1);
w=y(:,2);
plot(z,w) % the 1st component of solution 
xlabel('y_1')
ylabel('y_2') 
pause
end

end

t=t';
d=[t,y];
save 'data' d /ascii;

return