Recursive algorithms are derived to compute the generalized frequency response function matrix of multi-input multi-output (MIMO) non-linear systems as an analytical map from both non-linear differential equation models and NARX (Non-linear Auto Regressive Models with eXogenous inputs) models of the system. The algorithm is computationally compact and exposes the explicit relationship between the model parameters and the elements of the generalized frequency response function matrix and can thus provide important insights into the behaviour of non-linear systems.