magnetically dispersive property. Therefore, in order to have a classical treatment of electrodynamics in a magnetodielectric medium, we start with a Lagrangian for the total system (medium + electromagnetic field ) which is the sum of three parts $$\label{c1}
 L(t)= L_{res}+ L_{em}+L_{int}$$ where $L_{res}$ is the part related to the degrees of freedom of the medium and is the sum of the Lagrangians of the $E$ and $M$ fields $$\label{c2}
 L_{res}= L_{e}+L_{m}$$ where $$\label{c2.1}
L_{e}=\int _0^\infty d\omega \int d^3r \left[ \frac{1}{2}
 \dot{\vec{X}}_\omega \cdot  \dot{\vec{X}}_\omega
 -\frac{1}{2}\omega^2
 \vec{X}_\omega \cdot  \vec{X}_\omega\right]$$ and $$\label{c2.2}
L_{m}=\int _0^\infty d\omega \int d^3r \left[ \frac{1}{2}
 \dot{\vec{Y}}_\omega \cdot  \dot{\vec{Y}}_\omega
 -\frac{1}{2}\omega^2
 \vec{Y}_\omega \cdot  \vec{Y}_\omega\right]$$ Here the fields $ \vec{X}_\omega $ and $\vec{Y}_\omega $ are the dynamical variables of the $E$ and $M$ fields, respectively .

In (\[c1\]) $L_{em}$ is the contribution of the electromagnetic field in the Lagrangian of the total system $$\label{c3}
L_{em}= \int d^3r \left[ \frac{1}{2}\varepsilon _0
\vec{E}^2-\frac{\vec{B}^2}{2\mu_0}\right]$$ and $L_{int}$ is the part describing the interaction of the electromagnetic field with the medium $$\begin{aligned}
\label{c4}
&&L_{int}= \int _0^\infty d \omega \int d^3 r \int d^3 r' f_{ij}(
\omega , \vec{r} , \vec{r'}) E^i(\vec{r},t) \vec{X}_\omega
^j(\vec{r'} ,t)+\nonumber\\
&&\int _0^\infty d \omega \int d^3 r \int d^3 r' g_{ij}( \omega 