Years ago, I took a course related to this. The textbook was written by E. J. Mastascusa, who seems to be alive and well with his own web site.
Anyway, as I recall, we went over solutions of linear problems - sources, R, L, C, transformers, and dependent sources. All of these could be built out of standardized "little" matrices.
From a description file, the parts are given values and associated nodes. The values are populated into the part's matrix, and the matrices are superimposed into one large matrix according to their nodes. Independent sources are attached in a vector - again, according to their nodes.
Invert the large matrix, multiply it by the source vector and viola', you have the node voltages. Of course, this is only good for linear systems. It's great for DC operating points and simple AC analysis (Just put in the complex numbers). There was some way to handle transients using Rung-Kutta, but I've forgotten how.
I'm pretty sure the guys at Berkly came up with an auto assembling matrix as well, but I don't know how they would handle the nonlinear equations. Some very simple systems, like a battery, resistor, and diode, cannot be solved in closed form, so I think they must use an iterative solution.