function [tout,xout] = rk4(fun,t,x0)
% Function to demonstrate the 4th order Runge-Kutta Method for solution of IVP
% fun: is a function handle describing a first order ODE of the form
%       dx/dt = f(x,t)
% t: is a vector of time points. Assuming even spacing h = t(2)-t(1).
% x0: is the initial condition x(t=0)

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

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) + h*(k1+2*k2+2*k3+k4)/6;
end
