This method is taken from Taylor's Classical Mechanics, in the "Two-Body Central-Force Problems" section. This goes way more in-depth than it needs to, and assumes uniform densities of the stars for ease of calculation (though this is not necessary).
tl:dr: The Lagrangian is independent of each star's angle off the center of mass, and independent of each star's rotational angle, so the total orbital and each star's rotational angular momentum are all independently conserved. This is achieved by the changing mass accounting for the change in rotational speeds. Note that this does not require tidal locking or circular orbits.
Take $\vec{R}$ to be the center of mass position.
$$
\vec{R}=\frac{m_1\dot{}\vec{r}_1+m_2\dot{}\vec{r}_1}{m_1+m_2}=\frac{m_1\dot{}\vec{r}_1+m_2\dot{}\vec{r}_2}{M}
$$
where $M≡m_1+m_2$.
We can subsequently define $\vec{r}_1=\vec{R}+\frac{m_2}{M}\vec{r}$ and $\vec{r}_2=\vec{R}-\frac{m_1}{M}\vec{r}$, where $\vec{r}=\vec{r}_1-\vec{r}_2$
The kinetic energy is
$$
T=\frac{1}{2}(m_1\dot{\vec{r}}^2_1+m_2\dot{\vec{r}}^2_2+\frac{2}{5}m_1 s_1^2\omega_1^2+\frac{2}{5}m_2 s_2^2\omega_2^2)
$$
where $s_i$ are the radii of the stars, and $\omega _i$ are the rotational velocities (not orbital). We can reduce this to
$$
T=\frac{1}{2}(m_1\dot{\vec{r}}^2_1+m_2\dot{\vec{r}}^2_2+\frac{2}{5}m_1 s_1^2\omega_1^2+\frac{2}{5}m_2 s_2^2\omega_2^2)
$$
$$
T=\frac{1}{2}[m_1(\dot{\vec{R}}+\frac{m_2}{M}\dot{\vec{r}})^2+m_2(\dot{\vec{R}}-\frac{m_1}{M}\dot{\vec{r}})^2]+\frac{1}{5}(m_1 s_1^2\omega_1^2+m_2 s_2^2\omega_2^2)
$$
$$
T=\frac{1}{2}[M\dot{\vec{R}}^2+\frac{m_1 m_2}{M}\dot{\vec{r}}^2]+\frac{1}{5}(m_1 s_1^2\omega_1^2+m_2 s_2^2\omega_2^2)
$$
This lets us define a new quantity, the reduced mass: $\mu≡\frac{m_1 m_2}{M}$. We finally get
$$
T=\frac{1}{2}M\dot{\vec{R}}^2+\frac{1}{2}\mu\dot{\vec{r}}^2+\frac{1}{5}(m_1 s_1^2\omega_1^2+m_2 s_2^2\omega_2^2)
$$
For the total Lagrangian, taking a potential energy $U=U(r)$, then, we obtain
$$
L=T-U=\frac{1}{2}M\dot{\vec{R}}^2+\frac{1}{2}\mu\dot{\vec{r}}^2+\frac{1}{5}(m_1 s_1^2\omega_1^2+m_2 s_2^2\omega_2^2)-U(r).
$$
We can see that since the Lagrangian is independent of $\vec{R}$, that $M\ddot{\vec{R}}=0$ or $\dot{\vec{R}}=const.$. This tells us total momentum is conserved, our first conservation law. This is because, in the closed system, $\dot{m_1}=-\dot{m_2}$, so $\dot{M}=0$.
Since $\dot{\vec{R}}=const.$, we can move into the CM rest frame, so $\dot{\vec{R}}=0$.
$$
L=\frac{1}{2}\mu\dot{\vec{r}}^2+\frac{1}{5}(m_1 s_1^2\omega_1^2+m_2 s_2^2\omega_2^2)-U(r).
$$
Let $\dot{\vec{r}}^2=\dot{r}^2+r^2\dot{\phi}^2$, $\omega_i ^2=\dot{\theta}_i^2$.
$$
L=\frac{1}{2}\mu(\dot{r}^2+r^2\dot{\phi}^2)+\frac{1}{5}(m_1 s_1^2\dot{\theta}_1^2+m_2 s_2^2\dot{\theta}_2^2)-U(r).
$$
The Lagrangian is independent of $\phi$, so we again obtain a conservation equation.
$$
\frac{d}{dt}\frac{\partial L}{\partial \dot{\phi}}=0
$$
$$
\frac{\partial L}{\partial \dot{\phi}}=\mu r^2\dot{\phi}=const=l_{orbit}
$$
This tells us that orbital angular momentum is conserved. Apparently
$$
\mu r^2 \ddot{\phi} + 2 \mu r \dot{r} \dot{\phi} + \dot{\mu} r^2 \dot{\phi} = 0.
$$
If you were curious, $\dot{\mu}=\frac{\dot{m}_2 (m_1-m_2)}{M}=\frac{\dot{m}_1 (m_2-m_1)}{M}$.
Again, since the Lagrangian is independent of both $\theta _1$ and $\theta _2$, each star's individual rotational angular momentum is conserved.
We can find that
$$
\frac{2}{5}m_i s_i ^2 \dot{\theta}_i = const = l_{rot,i}
$$
and apparently
$$
m_i s_i ^2 \ddot{\theta}_i + 2 m_i s_i \dot{s_i} \dot{\theta}_i + \dot{m}_i s_i ^2 \dot{\theta}_i = 0,
$$
giving us our last constraint equation. Solving the Lagrangian for $r$ tells us the equations of motion.
There is another nice way of seeing this mathematically. It is not too hard to show that in the body frame, there are two conserved quantities: the square of the angular momentum vector
$$
L^2 = L_1^2 + L_2^2 + L_3^2
$$
and the rotational kinetic energy, which works out to be
$$
T = \frac{1}{2}\left( \frac{L_1^2}{I_1} + \frac{L_2^2}{I_2} + \frac{L_3^2}{I_3} \right).
$$
(Note that the angular momentum $\vec{L}$ itself is not conserved in the body frame; but its square does happen to be a constant.)
We can then ask the question: For given values of $L^2$ and $T$, what are the allowed values of $\vec{L}$? It is easy to see that $L^2$ constraint means that $\vec{L}$ must lie on the surface of a sphere; and it is almost as easy to see that the $T$ constraint means that $\vec{L}$ must also lie on the surface of a given ellipsoid, with principal axes $\sqrt{2TI_1} > \sqrt{2T I_2} > \sqrt{2T I_3}$. Thus, the allowed values of $\vec{L}$ must lie on the intersection of a sphere and an ellipsoid. If we hold $L^2$ fixed and generate a bunch of these curves for various values of $T$, they look like this:
Note that for a given value of $L^2$, an object will have its highest possible kinetic energy when rotating around the axis with the lowest moment of inertia, and vice versa.
Suppose, then, that an object is rotating around the axis of its highest moment of inertia. If we perturb this object so that we change its energy slightly (assuming for the sake of argument that $L^2$ remains constant), we see that the vector $\vec{L}$ will now lie on a relatively small curve near its original location. Similarly, if the object is rotating around its axis of lowest inertia, $\vec{L}$ will stay relatively close to its original value when perturbed.
However, the situation is markedly different when the object is rotating about the intermediate axis initially (the third red point in the diagram above, on the "front side" of the sphere. The contours of slightly perturbed $T$ near this point do not stay near the intermediate axis; they wander all over the sphere. There is therefore nothing keeping $\vec{L}$ from wandering all over this sphere if we perturb the object slightly away from rotating about this axis; which implies that an object rotating about its intermediate axis is unstable.
Best Answer
The angular momentum is conserved by this equation because it is derived from $$ \frac{\rm d}{{\rm d}t} \mathbf{L}_C = \sum \boldsymbol{\tau} $$
See Derivation of Euler's equations for rigid body rotation post for details. The derivative of angular momentum is zero when the torques are zero and thus $\mathbf{L}_C$ is constant.
I think the most likely scenario is that the numeric method will not conserve total energy, neither will it conserve angular momentum. There are some integration methods (called symplectic) that do conserve those quantities and Runge-Kutta is not one of them.
PS. NASA published an analytic solution to the above problem for some special cases (see this pdf report) and those analytic solutions do conserve angular momentum.