The simplest way to model nonlinear behavour is in the time domain.
One way to do it in the frequency domain if the response is for periodic (but non-sinusoidal, and amplitude dependent) is the Harmonic Balance Method (HBM or IHBM).
On the other hand, some nonlinear systems can generate responses that are unrelated to anything in a "linearised" (small-amplitude) model of the system or the frequency of the input force, and trying to model that in the frequency domain is a seriously hard problem.
The output from any frequency domain model tends to be restricted by the assumptions you made in setting up the model, and if those don't correspond to what the real system can do the "answers" can range from useful approximations to completely wrong. You need some independent information (either a time domain model, or some measurements) to validate what you get.
There is no particular problem in having a large (even infinite) frequency range in a FD model, but (as your OP implied) you need to careful about aliasing etc when converting between the time and frequency domains. This is often an artefact of using sampled data to represent the time domain, not something that is intrinsic to "real world" analog nonlinear systems.