The non-linear output frequency response functions (NOFRF), as an extension of the linear frequency response function (FRF) in the non-linear case, has been applied to weakly non-linear system study and engineering structural health monitoring (SHM). The computation of NOFRFs requires first solving a series of linear ordinary difference equations, i.e., generalized associated linear equations (GALEs), and then obtaining the system's results of each order according to the definition of NOFRFs in the frequency domain. However, in practical applications, the solution of GALEs often requires the aid of numerical integration. Therefore, accurate numerical computation of GALE is the first task in system analysis using NOFRFs. In our study, two different numerical methods are proposed for solving the system of linear differential equations of GALEs. The first computational method involves solving the GALEs of each order using a Recursive Computational Method (RCM). The second approach transforms the problem of solving GALEs into state-space equations, which are then solved using the integral solver of numerical computation software (e.g., MATLAB). This method is referred to as the coupled computational method (CCM). Finally, we compare the results of the two methods for computing NOFRFs using a non-linear differential equation (NDE) model with a fourth-order nonlinear term as an example. The final results show that the two methods give consistent results for low order NOFRFs. However, for higher order NOFRFs, CCM produces more accurate results than RCM. This provides ideas for calculating NOFRFs by GALE in nonlinear systems and also provides an important theoretical basis for calculating NOFRFs in multiple-input multiple-output (MIMO) systems.