Netlist: Math helps ... In the end the solution of the analog subsystem boils down to

(G - D) * V = I

with G being the conductance matrix, D a diagonal matrix with total conductance on the diagonal elements, V the net voltage vector and I the current vector.

By using solely two terminal devices, we can simplify the whole calculation significantly. A BJT now is a four terminal device with two terminals being connected internally.

The system is solved using an iterative approach:

G * V - D * V = I 

assuming V=Vn=Vo

Vn = D-1 * (I - G * Vo)

The above was already used in the code. However, with a pure two terminal device approach going forward we can e.g. further optimize. G typically is a sparse matrix.