In this article, we investigate how the entrywise application of a non-linear function to symmetric orthogonally invariant random matrix ensembles alters the spectral distribution. We treat also the multivariate case where we apply multivariate functions to entries of several orthogonally invariant matrices; where even correlations between the matrices are allowed. We find that in all those cases a Gaussian equivalence principle holds, that is, the asymptotic effect of the non-linear function is the same as taking a linear combination of the involved matrices and an additional independent GOE. The ReLU-function in the case of one matrix and the max-function in the case of two matrices provide illustrative examples.