How to model a circuit with resonant inductive coupling? I need to model (write the equations for) a circuit with resonant inductive coupling, such as the one in the picture below.

The picture is taken from the Wikipedia article on Resonant inductive coupling.
I need to write the differential equations that define this system and thus to calculate how much current flows through $Rl$, given $Vs$.
My problem is that I have zero clue where to start. I've only solved basic circuits so far (mainly DC). Even if I manage to calculate an equation for the left image (and perhaps solve it numerically), how can I use that for the right one?
My thoughts so far:
Using KVL for the first loop in the left circuit
$$V_s=V_{C_S}$$
And then applying this to the second loop:
$$ V_{C_S}=V_{R_S}+V_{L_S} $$
Then I substitued the formulas I knew for each component:
$$\frac{d^2 q}{dt^2}L_S+\frac{d q}{dt}R_S = \frac{q}{C_S}$$
Is my way of thinking good? Or did I neglect something that shouldn't have been?
And if It's good, after plugging this into a numerical calculator, how will I be able to tell the current in $Rl$?
Any given help would be greatly appreciated!
 A: I think you are on the right track. Just, no need to solve it numerically: second order linear equations with constant coefficients have a simple analytical solution.
The circuit on the left and the one on the right are variants of the RLC circuit, which is a damped oscillator. Try to solve the equation you have written with the condition that ${q \over Cs} = Vs = V_0\sin(\omega t)$. You will see that, depending on $\omega$, the circuit on the right responds in different ways.
The link between the two circuits is in the inductors. If the two circuits are close enough, the magnetic field produced by the inductor on the right will partly reach the inductor on the left. When the magnetic field changes, it will produce a voltage difference on the inductor on the left.
You can model this by writing that the voltage difference at $L_r$ is $V_{L_r} = -k L_s {dI_{L_s} \over dt}$, where $k$ is a coupling constant that depends on the distance between the two coils.
The above is true for any such circuit. If the circuit is also resonant, I believe this means that the resonant frequencies of the circuit on the right and on the left must be the same. But here ends my knowledge on the topic.

Addendum
Reading the comments, I see that an explanation on how to solve simple AC circuits in general might be required.
The differential equation describing an AC circuit is usually a second (or first) order linear ODE with an oscillating forcing term. It often has the form
$$a {d^2 y \over dt^2} + b {dy \over dt} + c y = A_0 cos(\omega t)$$
The oscillating term comes of course from the oscillating current/voltage of the alternating generator, and $\omega$ is the frequency of the generator.
This equation is very common in physics. Any time you talk about oscillations or resonance, chances are that this equation is behind it. The good news is that it is very easy to solve. You can look for sums of $sin$ and $cos$ functions, but if you are a bit familiar with complex numbers, they can save the day.
We can write the oscillating term as $A_0 \cos(\omega t) = \Re(A_0 e^{i\omega t})$.
Look for a solution in the form
$$y(t) =\Re(A e^{i(\omega t + \phi)})$$
where $A$ and $\phi$ are two real numbers to be determined. From now on I will not write "$\Re$" anymore and work just with complex numbers. You just need to remember to take the real part of the result after you have finished the calculation (you can check that this is indeed appropriate). Differentiating leads to
$${dy \over dt} = i\omega A e^{i(\omega t + \phi)} = i\omega y(t)$$
and
$${d^2y \over dt^2} = -\omega^2 y(t)$$
Substituting this in the differential equation leads to
$$-a\omega^2 y(t) + ib\omega y(t) + cy(t) = A_0 e^{i\omega t} = {A_0 \over A} e^{-i\phi} y(t)$$
$y(t)$ cancels in all the terms and you are left with an algebraic equation
$$A e^{i\phi} = \frac{A_0}{-a \omega^2 + ib\omega + c}$$
Defining $\omega_0 = \sqrt{c/a}$, the -so called- resonance frequency, this becomes
$$A e^{i\phi} = \frac{A_0/a}{(\omega_0^2 - \omega^2) + i(b/a)\omega}$$
Now, if you take the absolute value of both sides, you get
$$A = \frac{A_0/a}{\sqrt{(\omega_0^2-\omega^2)^2 + (b/a)^2 \omega^2}}$$
while the complex phase is
$$\phi = \tan^{-1}\left( \frac{(b/a)\omega}{\omega^2-\omega_0^2}\right)$$
And the final result is
$$y(t) = A \cos(\omega t + \phi)$$
More information and discussion on damped oscillators, for example here
