How do Circuit Simulators actually work?
I've examined the code of the Falstad simulator in some detail.  For circuits which consist only of linear components like resistors, switches, and voltage sources (things like logic-gate outputs are considered ground-connected voltage sources for purposes of the simulation) the simulator regards each circuit node, voltage source (connecting two nodes), or wire (likewise) as defining a linear equation and a variable, such that the number of equations and number of variables are always equal.  For a circuit node, the variable is the voltage of the node, and the equation computes the total current flowing through it equal to the total current injected by any current sources.  For a voltage source or wire (a wire being handled as a voltage source where the potential difference is zero), the equation sets the voltage difference between the two circuit node voltages equal to the required voltage difference, and the variable is the amount of current flowing through the voltage source from one node to the other.

Things like current sources and resistors are not associated with resistors or variables.  Instead, current sources increase the total current required for one circuit node (remember each circuit node has an equation which evaluates the total current flowing in and out) and decrease it for the other.  Resistors are a little trickier: for each endpoint's equation, the resistor adds terms for the node voltage of each endpoint.

A 100-ohm resistor connecting nodes 1 and 2, for example, would say that each volt increase on node 1 will decrease the current flowing into node 1 by 0.01 amps and increase the current flowing into node 2 by a like amount.  Likewise, each volt increase on node 2 would increase the current flowing into node 1 by 0.01 amps and decrease the current flowing into node 2 by a like amount.

Consider a circuit with a 10 volt supply connecting nodes 1 and 5, and 100 ohm resistors connecting node 1 and 2, 2 and 3, 2 and 4, and 3 and 4.  Assume further that there's a ground icon on node 1.  Thus:

neg ---+-1---R100---2---R100---3---100---4---pos
      gnd           |                    |
                    +---------100--------+


There would be two "voltage sources": the ground lead and the 10 volt supply (which are regarded as equation/variable 5 and 6, respectively).  The equations would thus be:

-X1*0.01                             +X5 -X6 = 0  Node 1
+X1*0.01 -X2*0.01          +X4*0.01          = 0  Node 2
         +X2*0.01 -X3*0.01 +X4*0.01          = 0  Node 3
         +X2*0.01          -X4*0.01      +X6 = 0  Node 4
-X1*1                                        = 0  Volts 5 (voltage between 1 and gnd)
-X1*1                      +X4*1             = 10 Volts 6 (voltage between 1 and 4)


This system of equation may be represented as an NxN matrix plus an N item array.  Each equation is represented by a row in the matrix, with values on each row representing the coefficients of each variable.  The right-hand side of each equation is stored in the separate array.  Before solving the equations, one will know the net current flowing into each node (zero in this case), and the voltage difference between pairs of nodes connected by voltage sources.  Solving the equations will yield the voltage at each node and the current flowing through each voltage source.

If the circuit contains capacitors, each of those will be regarded as a voltage source in series with a low-value resistor; after each simulation step, the voltage source will be adjusted according to the amount of current that flowed through it.  Inductors will be regarded as high-value resistors which feed current into one and and take it out the other (the amount of current being adjusted according to the voltage across the resistance).  For both capacitors and inductors, the value of the resistance will be controlled by the amount of time represented by a simulation step.

More complex circuit elements like transistors are regarded as combinations of voltage sources, current sources, and resistors.  Unlike the simpler circuit elements which let everything get processed once per simulation time step, elements like transistors compute their effective resistances etc. based upon the voltages and currents they're seeing, evaluate all the resulting equations, and re-evaluate their resistance based upon the new voltages and currents, re-evaluate the equations, etc. in an effort to reach an equilibrium where their effective resistance is as it should be for the voltage and current the transitor is seeing.

The Falstad simulator can be decently fast for moderate-sized circuits which consist entirely of "linear" elements.  The time to repeatedly solve a system of equations is pretty reasonable if the only thing that changes are the right-side coefficients.  The time gets much slower if the left side changes (e.g. because a transistor's effective resistance goes up or down) because the system has to "refactor" the equations.  Having to do refactor the equations multiple times per simulation step (may be necessary with transistors) makes things slower yet.

Using one big matrix for everything is not a good approach for large simulations; even though the matrix will be fairly sparse, it will take up space proportional to the square of the number of nodes plus voltage sources.  The time required to solve the matrix on each simulation step will be proportional to the square of matrix size if refactoring is not required, or to the cube of matrix size if refactoring is required.  Nonetheless, the approach does have a certain elegance when it comes to showing the relationship between a circuit and a system of linear equations.