I am writing a program where the user input $n$ variables, their initial values and differential equations. They may be non-linear.
To find the value of the variables at a time $(t_0 + T)$, I use 4th order Runge-Kutta method with fixed time step $h$.
What is the order of the error (both local and global) the user should expect?
By referencing Wolfram, I say the order of the local error is $O(h^5)$. Which makes the global error $O(h^4)$. All in a univariable context.
In a multivariable context the local and global error will be, I suspect, $O(h^{5n})$ and $O(h^{4n})$, respectively. But, I suspect again, this is true only in a linear context.
What should be the expected error of the 4th order Runge-Kutta integration in a multivariable (non-linear) context?
No, the error is always $O(h^4)$ independent of the dimension of the problem. The constant in that asymptotic bound of course depends on problem and thus dimension.
As a first estimate, use $Lh=\sqrt[5]\mu=10^{-3}$ where $L$ is a Lipschitz constant (estimate) for the problem.