Start off ignoring gravity. The spin axis is horizontal? Well then, you have an L vector. Very simple, nothing happens.
Now turn on gravity. This pulls down on the centre of mass, which is elsewhere from the pivot. That gives you an N vector - cross product of force with location relative to the pivot. This is perpendicular to the axis - so perpendicular to L. This is the change in L, according to N=dL/dt. In an arbitrarily small change in time, dt, we find dL will not change the length of L but will change its directions. Repeat indefinitely for every dt in a finite time interval t_1 to t_2.
I did not follow your derivation of $\frac{{\rm d}\mathbf{I}}{{\rm d}t}$. In most textbooks it is evaluated as follows $$\frac{{\rm d}\mathbf{I}}{{\rm d}t} =\boldsymbol{ \omega } \times \mathbf{I} = \begin{vmatrix} 0 & -\omega_z & \omega_y \\ \omega_z & 0 & -\omega_x \\ -\omega_y & \omega_x & 0 \end{vmatrix} \begin{vmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{vmatrix} $$
with the added caveat that $\mathbf{I}$ depends on the orientation of the body. The orientation may be tracked using Euler angles, Quaternions or just the 3×3 rotation matrix $\mathbf{R}$. Either way the end result is that the mass moment of inertia tensor needs to be calculated at every instant from the MMOI in the body frame
$$ \mathbf{I} = \mathbf{R}\,\mathbf{I}_{\rm body} \,\mathbf{R}^\top $$
In the end you have the equations of motion
$$ \left. \boldsymbol{\tau} = \mathbf{I}\, \boldsymbol{\dot{\omega}} + \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega}\;\; \right\} \;\; \boldsymbol{\dot{\omega}} = \mathbf{I}^{-1}\left(\boldsymbol{\tau} - \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega} \right) $$
It is also common to express the above in terms of angular momentum in the following algorithm. Each integration step is given the rotation matrix $\mathbf{R}$ and momentum vector $\boldsymbol{L}$
$$
\begin{array}{c|cc}
\text{Step} & \text{Calculation} & \text{Notes}\\
\hline 0 & \mathbf{I}=\mathbf{R}\mathbf{I}_{{\rm body}}\mathbf{R}^{\top} & \text{MMOI in world coorinates}\\
1 & \boldsymbol{\omega}=\mathbf{I}^{-1}\boldsymbol{L} & \text{Extract rotational vector}\\
2 & \dot{\mathbf{R}}=\boldsymbol{\omega}\times\mathbf{R} & \text{Change in rotation}^\star\\
3 & \dot{\boldsymbol{L}}=\boldsymbol{\tau}(t,\mathbf{R},\boldsymbol{\omega}) & \text{Change in momentum due to torque }\boldsymbol{\tau}
\end{array}
$$
* Note: When integrating the rotation matrix $\mathbf{R}$ using Runge-Kutta the result of $\mathbf{R} \rightarrow \mathbf{R} + h \dot{\mathbf{R}}$ is no longer a rotation matrix and the solution will quickly degrade in accuracy.
So instead, people often use quaternions $\boldsymbol{\hat{q}} = \pmatrix{ \boldsymbol{q}_{\rm v} & q_{\rm s}} $ that describe the rotation as $$ \mathbf{R} = \mathbf{1} + 2 q_{\rm s} [ \boldsymbol{q}_{\rm v}\times] + 2 [ \boldsymbol{q}_{\rm v} \times][ \boldsymbol{q}_{\rm v} \times] $$ where $[ \boldsymbol{q}_{\rm v} \times] = \begin{vmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{vmatrix}$ is the 3×3 cross product matrix operator of the vector part of the quaternion $\boldsymbol{q}_{\rm v}$.
The derivative of the quaternion is defined as $$ \dot{\boldsymbol{\hat{q}}} = \frac{1}{2} \pmatrix{ -\boldsymbol{\omega}^\top \boldsymbol{q}_{\rm v} \\ q_{\rm s} \boldsymbol{\omega} + \boldsymbol{\omega} \times \boldsymbol{q}_{\rm v} }$$
But often still people get this step wrong, because integrating the above $\boldsymbol{\hat{q}} \rightarrow \boldsymbol{\hat{q}} + h \dot{\boldsymbol{\hat{q}}}$ still brakes the rotation representation.
The proper way to take an integration step with quaternions is as follows. Given $\boldsymbol{\hat{q}} = \pmatrix{\boldsymbol{q}_{\rm v} & q_{\rm s}}$ and $\boldsymbol{\omega}$ vector
$$ \pmatrix{\boldsymbol{q}_{\rm v} \\ q_{\rm s}} \rightarrow \begin{vmatrix}
\cos(\tfrac{\theta}{2} ) & -\sin(\tfrac{\theta}{2} ) \boldsymbol{z}^\top \\
\sin(\tfrac{\theta}{2} ) \boldsymbol{z} & \cos(\tfrac{\theta}{2} ) + \sin(\tfrac{\theta}{2} ) [\boldsymbol{z}\times]
\end{vmatrix} \pmatrix{\boldsymbol{q}_{\rm v} \\ q_{\rm s}} $$
where $\theta = h \| \boldsymbol{\omega} \|$ is the step angle and $\boldsymbol{z} = \boldsymbol{\omega}/\|\boldsymbol{\omega}\|$ is the step rotation axis.
The resulting quaternion still represents rotations always and does not drift away like other formulations do.
Best Answer
$\vec\omega = I^{-1} \vec L$, and $\vec L$ is constant in the absence of external forces. The bit that I think you're missing is that $I$ rotates with the rigid body, so it is not constant in general and neither is $\vec\omega$.
I played with your online example, and the angular velocity does seem to always remain constant when I'm not poking the block, which is consistent with an inertia tensor that's a scalar multiple of the identity. I've never used Unity, but it seems to support arbitrary tensors of inertia via Rigidbody.inertiaTensor and Rigidbody.inertiaTensorRotation, so maybe you just need to set those properly.
If you end up having to roll your own physics, here's some naive, untested, and potentially wrong pseudocode: