Standard way to do Runge-Kutta (4th order) for coupled ODE's in Python?

Question

I am somewhat familiar with using RK4 for coupled ODE's, I found a very elegant way (in my opinion) to utilize it in Python, like so:

def RungeKutta4_2(x, y, z, dx, dydx, dzdx):
    
    k1 = dx*dydx(x, y, z)
    h1 = dx*dzdx(x, y, z)
    
    k2 = dx*dydx(x+dx/2., y+k1/2., z+h1/2.)
    h2 = dx*dzdx(x+dx/2., y+k1/2., z+h1/2.)
    k3 = dx*dydx(x+dx/2., y+k2/2., z+h2/2.)
    h3 = dx*dzdx(x+dx/2., y+k2/2., z+h2/2.)
    
    k4 = dx*dydx(x+dx, y+k3, z+h3)
    h4 = dx*dzdx(x+dx, y+k3, z+h3)

    y = y + 1./6.*(k1+2*k2+2*k3+k4)
    z = z + 1./6.*(h1+2*h2+2*h3+h4)
    x = x + dx
    
    return x, y, z

where x is the independent variable, and y and z are the dependent variables; I like this way of setting it up because the method itself automatically increments x by dx and the function just needs to be called in a while loop; the differential equations themselves are declared in their own functions. For example:

$$ \frac{dy}{dx} = z$$ $$ \frac{dz}{dx} = -y$$

def dydx(x,y,z):
    return z
def dzdx(x,y,z):
    return -y

while x <= x_end:
    x, y, z = RungeKutta4_2(x,y,z,dx,dydx,dzdx)
    x_values.append(x)
    y_values.append(y)
    z_values.append(z)

And then of course y_values and z_values could be plotted against x_values.

Is this the standard way to use RK4 for coupled ODE's? I am currently working on a paper I wish to publish, and my advisor mentioned that it is much more desirable to declare the ODE's in their own array, where dy/dx might be the first element and dz/dx is the second.

I'm having a lot of trouble figuring out why this would be "better", and if "my" way of using RK4 has any noticeable issues that my professor's way would fix?