function [ x_out, y_out ] = rk4_integrator(y0, xmin, xmax, N, f)
    dim = size(y0,2);
    y_out = zeros(N,dim);
    h = (xmax-xmin)/(N-1); 
    
    x = xmin;
    for i = 1:N
        y_out(i,:) = y0;
        y0 = rk4(y0,x,h,f);
        x = xmin+h*i;
    end
    
    x_out = linspace(xmin, xmax, N);
end

