Measure and simulate non-linear distortion
Nonlinear systems are very hard to classify and there's no unified theory like for linear systems. In general, you cannot measure/identify non-linear systems in finite time. There are  some specific classes of nonlinear systems which allow for identification or at least approximation. 

You already named a trivial one: The memoryless nonlinear system. In this case you can learn the nonlinear transfer function (i.e. your waveshaping function) by just looking at the map from input to output.

Any interesting real system does have memory however, and there we can start with a specific model for the system and adapt the parameters until the model output and the real output match. There are nonlinear optimization techniques that will help to find the proper set of parameters.

If you know little about the internals of the system you want to simulate, then you can use a generic model. The most popular such model is described by the Volterra series and requires that the system is expandable in a series similar to the taylor series but involving convolution terms and that this series converges quickly enough (and in the entire domain you're interested in) to be practically useful. The theory of Volterra kernels delivers the mathematical tool to find the series expansion from input/output observations without the need for non-deterministic optimization techniques.

So summing up, your best chance to really capture the behavior of a nonlinear system is to understand the system in detail and model its components. The black box approach is very difficult and only works under very specific conditions.