function [varargout] = rk4(x0,y0,z0,beta,sigma,ro,h,steps)
%This function is the numerical ode solver using the Runge Kutta 4th order
%method.

        %The output variables are initialized for speedy execution 
        x_fin = zeros(1,steps);
        y_fin = zeros(1,steps);
        z_fin = zeros(1,steps);
        
        x_fin(1) = x0;
        y_fin(1) = y0;
        z_fin(1) = z0;
        %The numerical integration starts from here
        for count = 2:steps
            f1 = lor_eq(x_fin(count-1),y_fin(count-1),z_fin(count-1),beta,sigma,ro);
            f2 = lor_eq(x_fin(count-1)+0.5*h*f1(1),y_fin(count-1)+0.5*h*f1(2),z_fin(count-1)+0.5*h*f1(3),beta,sigma,ro);
            f3 = lor_eq(x_fin(count-1)+0.5*h*f2(1),y_fin(count-1)+0.5*h*f2(2),z_fin(count-1)+0.5*h*f2(3),beta,sigma,ro);
            f4 = lor_eq(x_fin(count-1)+h*f3(1),y_fin(count-1)+h*f3(2),z_fin(count-1)+h*f3(3),beta,sigma,ro);
            vec = [x_fin(count-1);y_fin(count-1);z_fin(count-1)] + h*(f1 + 2*f2 + 2*f3 + f4)/6;
           
            
            x_fin(count) = vec(1);
            y_fin(count) = vec(2);
            z_fin(count) = vec(3);
        end
        varargout(1) = {x_fin};
        varargout(2) = {y_fin};
        varargout(3) = {z_fin};

end

        
