function [t,x] = runge(field, t0, T, x0, h)
%RUNGE Implements the Runge one-step ODE solver
%
% Parameters:
%  field  -- Right hand side function of ODE system: x'=f(t,x)
%  t0     -- Initial time
%  T      -- End time (T > t0)
%  x0     -- Initial value
%  h      -- Size of time step (h <= T-t0)
%
% Outputs:
%  t  -- [t0; t-0+h, t0+2*h; ...; t0+i*h; ...]
%  x  -- Column vector containing numerical solution at each time step

n = ceil((T-t0)/h);
t = t0+h*(0:n).';
x = ones(n+1,length(x0));
x(1,:) = x0;

for i=1:n
    k1 = feval(field, t(i), x(i,:).');
    k2 = feval(field, t(i)+0.5*h, x(i,:).'+0.5*h*k1);
    x(i+1,:) = x(i,:).' + h*k2;
end