function [x,y] = genericODE(x_1,x_f,h,y_1,deriv,odeMethod)
%[x,y] = genericODE(x_1,x_f,h,y_1,deriv,odeMethod);
% a generic computation scheme for integration of a system of ode (this 
% module eliminates the need to keep writing the same code over
% and over again, allowing for arbitrary derivatives and ode integration 
% methods)
%      x_1 --> initial independent variable
%      x_f --> final independent variable
%        h --> step size in independent variable
%      y_1 --> initial dependent variable(s) [column vector]
%    deriv --> function handle to system of ode
%odeMethod --> function handle to ode integration method
%        x <== vector of (x_1:h:x_f)' [note this is a column vector]
%        y <== solutions to ode(s) at x [each row is y(x)]
%
% Note that this code assumes column vector output for your deriv function,
% but you may have written it to give row vector output!  this is easy to
% remedy with your function handle deriv.  use the following scheme:
%
%  deriv gives column vectors:  deriv = @(x,y) foo(x,y,a,b,...)
%  deriv gives row vectors:     deriv = @(x,y) foo(x,y,a,b,...).'
% 
% odeMethod is function handle to the integration method, an example of
% which is given by:
%
% odeMethod = @(x,y,h,deriv) odeStep_euler(x,y,h,deriv);

% (c) 2008 Curt A. L. Szuberla
% University of Alaska Fairbanks, all rights reserved

% ver 1.00
% 3 February 2008

% initialize/pre-allocate independent and dependent variables
x = (x_1:h:x_f)';                  % each row is a step
y = zeros(length(x),length(y_1));  % each row is a solution y(x)

% initial conditions
y(1,:) = y_1';

% compute the ode solutions
for k = 2:length(x)
    % you should replace odeStep_euler, below, with the name of your own
    % Euler step function
    y(k,:) = odeMethod(x(k-1),y(k-1,:).',h,deriv);
end

return