function [tout,xout] = rk4sys(fun,t,x0)
% Function to demonstrate the 4th order Runge-Kutta Method for solution of
%  a system of IVP

n = length(x0);

h = t(2)-t(1);
tout = zeros(length(t),1);
tout(1)=t(1);

xout = zeros(length(t),n);
xout(1,:)=x0(:);

for i=1:length(t)-1
    k1 = fun(xout(i,:),tout(i))';
    k2 = fun(xout(i,:)+h*k1/2,tout(i)+h/2)';
    k3 = fun(xout(i,:)+h*k2/2,tout(i)+h/2)';
    k4 = fun(xout(i,:)+h*k3,tout(i)+h)';
    
    tout(i+1) = tout(i)+h;
    xout(i+1,:) = xout(i,:) + (k1+2.*k2+2.*k3+k4).*h./6;
end
