Apparent violation of energy conservation in an electromagnetic transformer system First of all I would like to say I am not a native English speaker, so there might be some grammar mistakes. Sorry in advance for that!
I have been a whole month trying to understand what is going on in a certain electromagnetic system, with no success. Maybe some of you can shed some light on it.
Picture the following system: There is an ideal 1:1 electric transformer. Both coils have inductance L, and mutual inductance takes the value $M=\sqrt{L_1 L_2}=\sqrt{L^2}=L$, as there are no losses and the whole magnetic flux created by each coil goes through the other one.
Now, on both sides of the transformer, we connect a resistor and a capacitor in the way shown in the picture:

The resistors have both resistance R, and the capacitor have both capacitance C.
The only difference between both sides is that the capacitor on the left hand side is initially charged with charge Q, while the one on the right hand side has no initial charge. Initially, there is no current going through the coils either.
As far as I know, the time evolution of the charge stored in the capacitors is given by the following two coupled differential equations, where $x_1$ is the charge in the left-hand side capacitor and $x_2$ is the charge in the right-hand side one, both as functions of time.
$$
L\frac{d^2x_1(t)}{dt^2}+M\frac{d^2x_2(t)}{dt^2}+R\frac{dx_1(t)}{dt}+\frac{x_1}{C}=0
$$
$$
L\frac{d^2x_2(t)}{dt^2}+M\frac{d^2x_1(t)}{dt^2}+R\frac{dx_2(t)}{dt}+\frac{x_2}{C}=0
$$
With initial conditions:
$$
\frac{dx_1}{dt}(0)=0 \ \ \ \ \ x_1(0)=Q
$$
$$
\frac{dx_2}{dt}(0)=0 \ \ \ \ \ x_2(0)=0
$$
I have checked this equations over and over and I honestly think they are right. However, the problem arises when I solve them. I have built a simulink model to do it, as well as a matlab code using ode45. When I run the script, I get the charge in the capacitors as functions of time. Taking the first derivatives, I get the intensity in both sides of the transformer.
Now, if I evaluate both the intensities and the stored charges in a big enough time from the start, they both tend to zero as expected. The energy initially stored in the left capacitor gets dissipated in the resistors within a few cycles. The problem is that if I calculate the dissipated energy, by doing the following integral:
$$
U_{diss}=\int_0^\infty R \left(\frac{dx_1(t)}{dt} \right)^2 dt + \int_0^\infty R \left(\frac{dx_2(t)}{dt} \right)^2 dt
$$
I get a bigger number than if I calculate the initially stored energy in the capacitor. This would mean that somehow the circuit dissipated more energy by Joule's effect than the initial energy stored in it. 
Obviously I am doing something wrong, but I have no idea about what could it be.
I am completely sure that I have solved the differential equations correctly, as I have checked the numerical method used many times. I have built several models using different software and I always get the same result, and a friend on mine that is an expert in numerical methods for differential equations have checked it too. So the problem must be in the differential equations themselves, not it their solution. But the more I check them, the more I am convinced they are right, so I have no idea about what to do.
To me, it is clear that the dissipated energy for large values of time should be the initially stored energy in the capacitors, never bigger.
Please help me with this. I really need to get it right.
Thanks in advance for your answers,
Alberto.
 A: There is nothing wrong with your equations that leads to non-conservation of energy. You can see this by deriving the energy conservation law directly from the equations of motion. For this, multiply the first equation by $dx_1/dt$, the second by $dx_2/dt$, and sum them. After collecting some terms into total derivatives, you obtain the following (dots are time derivatives):
\begin{equation}
\frac{d}{dt}\left( \frac{1}{2}L \dot{x_1}^2 + \frac{1}{2}L \dot{x_2}^2 + \frac{1}{2C} {x_1}^2 + \frac{1}{2C} {x_2}^2 + M \dot{x_1}\dot{x_2}  \right) = -R \dot{x_1}^2 - R \dot{x_2}^2
\end{equation}
You can identify what you'd like to call energy on the left-hand side. Integrating it and plugging in your initial conditions, you get
\begin{equation}
E(0) - E(\infty) = \frac{Q^2}{2 C} = \int \limits_0^\infty dt \hspace{3pt} R \dot{x_1}^2 + R\dot{x_2}^2 
\end{equation}
All this follows directly from the equations, without any reference to physics behind it. If the last equation is violated, it is a defect of either numerical integration scheme (less likely), or your ODE solver (more likely).
A: Most of the things you've done are correct except the evaluation of the last integral.
$\frac{d^2x}{dt^2}\neq(\frac{dx}{dt})^2$
So, when you put that as a function into your script, you will get wrong results! 
