function y = RK4(f,x,y0)
%
% RK4 uses the classical RK4 method to solve a first-order ODE 
% given in the form y' = f(x,y) subject to initial condition y0.
%
%  y = RK4(f,x,y0) where
%
%       f is an inline function representing f(x,y),
%       x is a vector representing the mesh points,
%       y0 is a scalar representing the initial value of y,  
%
%       y is the vector of solution estimates at the mesh points.

y = 0*x;    % Pre-allocate
y(1) = y0; h = x(2)-x(1);
for n = 1:length(x)-1,
    k1 = f(x(n),y(n));
    k2 = f(x(n)+h/2,y(n)+h*k1/2);
    k3 = f(x(n)+h/2,y(n)+h*k2/2);
    k4 = f(x(n)+h,y(n)+h*k3);
   y(n+1) = y(n)+h*(k1+2*k2+2*k3+k4)/6;
end