There is a growing interest in Virtual Analog modeling algorithms for musical audio processing designed in the Wave Digital (WD) domain. Such algorithms typically employ a discretization strategy based on the trapezoidal rule with fixed sampling step, though this is not the only option. In fact, alternative discretization strategies (possibly with an adaptive sampling step) can be quite advantageous, particularly when dealing with nonlinear systems characterized by stiff equations. In this article, we propose a unified approach for modeling capacitors and inductors in the WD domain using generic linear multi-step discretization methods with variable time-step size, and provide generalized adaptation conditions. We also show that the proposed approach for implementing dynamic (energystoring) elements in the WD domain is particularly suitable to be combined with a recently developed technique for efficiently solving a class of circuits with multiple one-port nonlinearities, called Scattering Iterative Method. Finally, as examples of application, we develop WD models for a Van Der Pol oscillator and a dynamic diode-based ring modulator, which use different discretization methods.