## 1 Introduction
PWL approximation of the nonlinear characteristics of circuit elements allows the problem to be solved by checking some regions (or hyper-regions) for solutions. However, as shown in numerical experiments [18], the PWL approximation often changes the number of solutions. Thus, to design a high-reliability circuit, it is necessary to find all solutions to the original nonlinear circuits. For this purpose, interval methods have been developed. The methods are capable of finding all DC solutions at a given accuracy [16; 19; 21; 23; 24]. However, the computational time of interval algorithms increases exponentially with the number of variables. Introducing various tests to interval algorithms (e.g., linear programming test) allows all DC solutions to be determined more efficiently [18; 19]. A similar approach is based on the idea of successive contraction, division, and elimination of certain hyper-rectangular regions in which solutions are sought. The essential point of this approach is the contraction method, which leads to obtaining bounds on the location of all solutions inside a hyper-rectangular region [20; 22]. Although all the methods are mathematically elegant, they are rather time-consuming and computationally intensive for large-scale circuits. Consequently, only relatively small or  
Figure 1: Schmitt trigger and the results of default DC analysis in IsSpice4 software. specific nonlinear circuits can be effectively analyzed. In addition, these methods usually require a specific hybrid description, apply to circuits belonging to a strictly defined class, and use simple models of semiconductor elements (e.g., the Ebers-Moll model of a bipolar transistor, the Ebers-Moll-type model (The Ebers-Moll model of bipolar transistors is composed of two diodes and two current-controlled current sources. The MOS transistors with a long channel can be characterized by a model having the same structure, called the Ebers-Moll-type model. This model is equivalent to the Shichman-Hodges model with the omission of some secondary phenomena [25; 26]) or the Shichman-Hodges model of an MOS transistor) [19; 22; 27]. The theoretical problem of finding all solutions to an electronic circuit of significant size is still unsolved and, according to many researchers, will remain unsolved and challenging for a long time to come [28]. Therefore, some methods focus on determining only a certain subset of these solutions. Many homotopy methods have been developed to calculate multiple solutions [5; 6; 11; 14; 27; 29; 30; 31; 32]. Most homotopy methods follow the multi-start or extended curve tracking (lambda-threading) approach. The multi-start approach selects multiple initial guesses at the homotopy parameter, lambda, equal to zero and follows the continuation paths to solutions at lambda = 1. The extended curve-tracing approach follows a single homotopy path past lambda = 1, hoping to reverse direction and pass through lambda = 1 multiple times. This approach may not find all solutions if the tracked path does not pass through all solutions, but it is efficient because it can avoid revisiting the same solution [29]. For example, [PERSON] [6] developed a SPICE-oriented Newton homotopy method based on the arc length approach, which chooses one of the node voltages in each positive feedback loop as a variable and traces the solution path using SPICE transient analysis. The work [27] brought a method using thread homotopy based on transistor states. The papers [11; 29] provide a comprehensive analysis of homotopy methods for calculating all DC solutions of nonlinear circuits and systems. Several sufficient conditions were derived to guarantee finding all solutions along a single homotopy path. An efficient method for finding operating points of nonlinear circuits containing neither voltage- nor current-controlled PWL resistors using the homotopy method and SPICE-oriented approach was proposed in [14]. Finding multiple DC operating points in memristor circuits using homotopy methods was studied in [31].