The simulation of nonlinear components is central to virtual analog simulation. In audio effects, circuit elements often include devices such as diodes and transistors, mostly operating in the strongly nonlinear regime. Mathematical models are of the form of systems of nonlinear ordinary differential equations (ODEs), and traditional integrators, such as the trapezoid and midpoint methods, can be employed as solvers. These methods are fully implicit, and require the solution of a nonlinear algebraic system at each time step, introducing further complications regarding the existence and uniqueness of the solution, as well as the choice of halting conditions for the iterative root finder. On the other hand, fast explicit methods such as Forward Euler, or explicit Runge-Kutta methods, are prone to unstable behaviour at standard audio sample rates, even at moderate amplitudes. For these reasons, in this work a family of linearlyimplicit schemes is presented. These schemes take the form of a perturbation expansion, making the construction of higherorder schemes possible. Compared with classic implicit designs, the proposed methods have the advantage of efficiency, since the update is computed in a single iteration, through the solution of a linear system of equations. Furthermore, the existence and uniqueness of the update are proven by simple inspection of the update matrix. Compared to classic explicit designs, the proposed schemes display stable behaviour at standard audio sample rates. In the case of a single scalar ODE, sufficient conditions for numerical stability can be derived, imposing constraints on the choice of the sampling rate. Several theoretical results are provided, as well as numerical examples for typical stiff equations used in virtual analog modelling.