What you call n coupled equations is in Jamal's reference a single equation

dy(x)/dx = f(x,y)

for a n-component quantity y=(y_1,....y_n).

The solution algorithm rk4 (p. 712) depends on n simply as a parameter. Whether this parameter has the value 3 or 100 makes no difference on the level of coding. Execution time should be propotional to n.

See numerical recipes page 710-716 (Numerical Recipes in C

The Art of Scientific Computing

Second Edition, CAMBRIDGE UNIVERSITY PRESS)

What you call n coupled equations is in Jamal's reference a single equation

dy(x)/dx = f(x,y)

for a n-component quantity y=(y_1,....y_n).

The solution algorithm rk4 (p. 712) depends on n simply as a parameter. Whether this parameter has the value 3 or 100 makes no difference on the level of coding. Execution time should be propotional to n.

One can use the RK4 for the system of ODEs the number of equations is not a question. I have written a code in MATLAB (same as Mathews Book) and only thing required in the code is the initial conditions. If you have 100 initial conditions then it will solve without explicitly writing ks

You just use the same integration algorithms, using linear combinations of vectors as the arguments to the function f in the equation y' = f'(y,t). There is no essential difference. With due care the same algorithms can be used in the case if differential in difference equations.