function xout = step_it(func, xin, par, stepsize, dim) % The function stepit is used to integrate a system of autonomous % differential equations. % % It performs a numerical integration over a time 'stepsize' of a % specified function, 'func', using a 4th order Runge-Kutta scheme % (see, for example, Press et al, 1986. Numerical Recipes in FORTRAN, % Cambridge University Press). % % Useage; % stepit('func', xin, par, stepsize, dim) % % where func is the name of the function to be integrated, xin is the % current state of the system, par contains any necessary parameters, % stepsize is the integration step size, and dim is the number of ODEs. % % You would normally include this function in a loop. x = xin; f = feval(func, x, par, dim); c1 = stepsize .* f; x = xin + c1 /2; f = feval(func, x, par, dim); c2 = stepsize .* f; x = xin + c2 /2; f = feval(func, x, par, dim); c3 = stepsize .* f; x = xin + c3; f = feval(func, x, par, dim); c4 = stepsize .* f; xout = xin + (c1 + 2.*c2 + 2.*c3 + c4)./6;