I have always taken for granted circuit simulators and I didn't spend time understanding now. I wish to understand better now.

Normally when you forward simulate ODE systems there is a single dynamic equation which you forward integrate(ex, simplest way is Eulers method). What I am wondering about is in the case of circuit design (SPICE, LT-SPICE, etc) how do you maintain the core constraints of the system? EX KCL and KVL (Kirchoff voltage/current law) constraints always need to be maintained and it's not clear how these constraints can be factored in when integrating system dynamics.