You might be interested in the free education edition of PERMAS. The Harmonic Balance Method (HBM) is available to compute nonlinear frequency response functions. See https://www.intes.de/edu
Some academic HBM examples are available in the permas4edu project on github. See
Years ego, I wrote a program to solve a second order differntial equation using Runge Kuta numerical method in BASIC language to simulate a spring-damper-mass system.
This method can be used for any nonlinear dynamic system as well as duffing equation that is a second order differntial equation and any scientific math software or programing language can be used to find the solution and visualize the results.
Here is a tutorial I wrote on how to obtain the frequency response using numerical continuation of periodic orbits: FrequencyResponse.pdf (utwente.nl). This uses MatCont which is a free Matlab continuation toolbox for bifurcation analysis.
A relatively simple approach involves using numerical methods (such as the 'ode45' function in MATLAB) to obtain the time history information of the Nonlinear Dynamic system, and then utilizing Fast Fourier Transform (FFT) in MATLAB to obtain the amplitude-frequency response, thereby obtaining the frequency response you need. Below is a simple document I've written along with the corresponding MATLAB program. You might consider trying it out to solve your own problem. If there are any incorrect considerations, I apologize and kindly ask for them to be pointed out.