function [rate] = dydt(k, y, t, dt ) 

% y & rate should be 1xD row vector, D is number of equations 
% returns "average" value of rate (dy/dt) over interval t to t + dt 
% given the function f_RHS from the original ODE dy/dt = f_RHS(y,t) 
% f_RHS must be available in a separate function file 

if k==1
% Euler. 

rate = f_RHS( y , t );

elseif k==2
% Modified Euler. 

dydt0 = f_RHS( y , t ); 
dydt1 = f_RHS( y + dt*dydt0 , t + dt ); 
rate = 0.5*(dydt0 + dydt1); 

elseif k==4
% 4th order Runge-Kutta. 

dydt0 = f_RHS( y , t ); 
dydt1 = f_RHS( y + 0.5*dt*dydt0 , t + 0.5*dt ); 
dydt2 = f_RHS( y + 0.5*dt*dydt1 , t + 0.5*dt ); 
dydt3 = f_RHS( y + dt*dydt2 , t + dt ); 
rate = dydt0/6 + dydt1/3 + dydt2/3 + dydt3/6; 

else 

rate = 0;

end