Numerical modeling of induced voltage

OK, I figured it out: Apply Faraday's induced voltage law on each of the 2 coils, using the magnetic flux due to each of the 2 coils, yielding 2*2=4 total contributions. And then solve the system of equations so the currents are consistent with the induced voltages.

I used my favorite optimizer/constraint solver, cvxpy to do it. Here's the Python code:

def solveForCurrents(i1_init, i2_init,  flux1dueto1, flux2dueto1, flux2dueto2, flux1dueto2, r1,r2,external_emf_1,external_emf_2, dt):
    """ Convention:
    flux2dueto1 is (coil2 number of turns) * (the surface integral of magnetic flux through coil 2, per amp through coil 1).
    external_emf_1 is a voltage in series with coil 1, e.g. the instantaneous voltage on a capacitor.
    r1 is the total resistance of coil 1, including equivalent series resistance (ESR) of capacitors etc.
    """
    import cvxpy
    flux1_init = flux1dueto1*i1_init + flux1dueto2*i2_init
    flux2_init = flux2dueto1 * i1_init + flux2dueto2 * i2_init
    i1 = cvxpy.Variable(1) #variables to be solved for in the system of equations
    i2 = cvxpy.Variable(1)
    flux1 = flux1dueto1 * i1 + flux1dueto2 * i2 #total flux through a coil is the sum of all sources. (The flux1dueto1 captures the effect of self-inductance)
    flux2 = flux2dueto1 * i1 + flux2dueto2 * i2
    induced_emf1 = -1 * (flux1 - flux1_init) / dt #Faraday's law: induced voltage = -1 * change in flux per time.
    induced_emf2 = -1 * (flux2 - flux2_init) / dt
    constraints = [(induced_emf1 + external_emf_1) == i1 * r1, #Ohm's law: V=IR
                   (induced_emf2 + external_emf_2) == i2 * r2]
    prob = cvxpy.Problem(cvxpy.Minimize(1), constraints=constraints) #nothing to minimize; just solving for the constraints
    prob.solve()
    return i1.value[0], i2.value[0]

#example usage. Would run this in a loop supplying i1 and i2 as the next timestep's i1_init and i2_init 
i1, i2 = solveForCurrents(i1_init=0, i2_init=0,  flux1dueto1=0.1, flux2dueto1=0.01, flux2dueto2=0.01, flux1dueto2=0.001, r1=1,r2=10,external_emf_1=1000,external_emf_2=0, dt=0.001)

In my case, to get the flux2dueto1 for a given coil geometry, I used the Biot-Savart law with a 1 amp current through coil 1, to find the magnetic field at a point on the face of coil 2. I repeated for a 20x20 grid across the face of coil 2, numerically integrating the results to get the integral of flux through the coil. (That surface-integral code isn't shown above; my implementation was pretty specific to my geometry.)

I can cache and reuse that same value of "flux per amp" to get the flux due to any value of current (as long as the geometry stays the same. If geometry changes then flux2dueto1 needs to be recalculated).