function [uout, tout, error ] = rk4(f_1, f_2, f_3, ti, tf, k, u0 )
dt = k;
nt_steps = (tf - ti)/dt;


vexact = @(t) -sin(2.*t) + t.^2 - 3;
%Initializing Arrays
t = linspace(0, 2, nt_steps);
x = zeros(1, nt_steps);
y = zeros(1, nt_steps);
z = zeros(1, nt_steps);


vbar = vexact(t);
x(1) = u0(1);
y(1) = u0(2);
z(1) = u0(3);

for i  = 1:(nt_steps - 1)
    k_1 = k*f_1(t(i),x(i),y(i),z(i));
    L_1 = k*f_2(t(i),x(i),y(i),z(i));
    m_1 = k*f_3(t(i),x(i),y(i),z(i));

    k_2 = k*f_1(t(i) + 0.5*dt, x(i) + ...
        0.5*k_1, y(i) + 0.5*L_1, z(i) + 0.5*m_1);
    L_2 = k*f_2(t(i) + 0.5*dt, x(i) + ...
        0.5*k_1, y(i) + 0.5*L_1, z(i) + 0.5*m_1);
    m_2 = k*f_3(t(i) + 0.5*dt, x(i) + ...
        0.5*k_1, y(i) + 0.5*L_1, z(i) + 0.5*m_1);


    k_3 = k*f_1(t(i) + 0.5*dt, x(i) + ...
        0.5*k_2, y(i) + 0.5*L_2, z(i) + 0.5*m_2);
    L_3 = k*f_2(t(i) + 0.5*dt, x(i) + ...
        0.5*k_2, y(i) + 0.5*L_2, z(i) + 0.5*m_2);
    m_3 = k*f_3(t(i) + 0.5*dt, x(i) + ...
        0.5*k_2, y(i) + 0.5*L_2, z(i) + 0.5*m_2);


    k_4 = k*f_1(t(i) + dt, x(i) + ...
        k_3, y(i) + L_3, z(i) + m_3);
    L_4 = k*f_2(t(i) + dt, x(i) + ...
        k_3, y(i) + L_3, z(i) + m_3);
    m_4 = k*f_3(t(i) + dt, x(i) + ...
        k_3, y(i) + L_3, z(i) + m_3);


    x(i+1) = x(i) + (k_1 + 2*k_2 + 2*k_3 + k_4)/6;
    y(i+1) = y(i) + (L_1 + 2*L_2 + 2*L_3 + L_4)/6;
    z(i+1) = z(i) + (m_1 + 2*m_2 + 2*m_3 + m_4)/6;
end
uout = [x;y;z];
tout = t;
error = max(max(abs(x(5) - vbar(5))));




end