No, Euler angles are a poor choice for representing the orientation. Apart from that not all orientations can be represented it complicates the calculations.
What you instead could do is to use quaternion representation of orientation or use a rotation matrix. That would make the maths more easy, the relation between angular velocity and orientation is $\omega q = {1\over 2}{d\over dt}\operatorname{vec}(q)$, where $\operatorname{vec}(q)$ is the vector part of the quaternion.
I'll focus my answer on computing the matrix exponential. Consider components of the unit vector $\hat{\omega}^a = \omega^a /|\omega| $. Splitting the exponential into odd/ even terms
$$ E=\sum_{n=0}^\infty \frac{\left( i |\omega|\hat{\omega}^a S^a \right)^n}{n!} = 1+\sum_{n=1}^\infty \frac{(-1)^n |\omega|^{2n}(\hat{\omega}^aS^a)^{2n}}{(2n)!}+i \sum_{n=0}^\infty \frac{(-1)^n |\omega|^{2n+1}(\hat{\omega}^aS^a)^{2n+1}}{(2n+1)!} $$
The components $jk$ of the $\operatorname{so}(3)$ rotation generator $a$ are $S^{ajk}=i\epsilon^{ajk} $. Consider
$$\left[ (\hat{\omega}^aS^a)^2\right]^{jl}=\hat{\omega}^a \hat{\omega}^b S^{ajk}S^{bkl}=\delta^{jl}-\hat{\omega}^j\hat{\omega}^l$$
The last equality follows$^\dagger $ by using the 'contracted epsilon identity' as well as the fact that $\hat{\omega}^a \hat{\omega}^a =1 $. We also note that RHS is a projection operator$^\ddagger$, hence
$$\left[ (\hat{\omega}^aS^a)^{2n}\right]^{jl}=\delta^{jl}-\hat{\omega}^j\hat{\omega}^l$$
$$\left[ (\hat{\omega}^aS^a)^{2n+1}\right]^{jl}=\left[ \hat{\omega}^aS^a \right]^{jl}$$
Substitute back into the sums
$$ E^{jl}=\delta^{jl}+\sum_{n=1}^\infty \frac{(-1)^n |\omega|^{2n}}{(2n)!}(\delta^{jl}-\hat{\omega}^j\hat{\omega}^l) +i \sum_{n=0}^\infty \frac{(-1)^n |\omega|^{2n+1}}{(2n+1)!}\left[ \hat{\omega}^aS^a \right]^{jl} $$
$$E^{jl}=\delta^{jl}+(\delta^{jl}-\hat{\omega}^j\hat{\omega}^l)\left(\cos(|\omega|)-1\right) +i\left[ \hat{\omega}^aS^a \right]^{jl}\sin(|\omega|)$$
$$E^{jl}=\hat{\omega}^j\hat{\omega}^l+(\delta^{jl}-\hat{\omega}^j\hat{\omega}^l)\cos(|\omega|) +i\left[ \hat{\omega}^aS^a \right]^{jl}\sin(|\omega|)$$
Those are the components of the matrix. In matrix form it's probably not going to be very pretty, so I'll leave it there.
$\dagger$ Like this
$$ \hat{\omega}^a \hat{\omega}^b S^{ajk}S^{bkl}=i^2\hat{\omega}^a \hat{\omega}^b \epsilon^{ajk}\epsilon^{bkl}=\hat{\omega}^a \hat{\omega}^b \epsilon^{ajk}\epsilon^{blk}=\hat{\omega}^a \hat{\omega}^b \left( \delta^{ab}\delta^{jl}-\delta^{al}\delta^{bj} \right) \\= \hat{\omega}^a \hat{\omega}^a \delta^{jl}-\hat{\omega}^j \hat{\omega}^l = \delta^{jl}-\hat{\omega}^j \hat{\omega}^l $$
$\ddagger$ By direct computation
$$ (\delta^{jl}-\hat{\omega}^j \hat{\omega}^l)(\delta^{lm}-\hat{\omega}^l \hat{\omega}^m)= \delta^{jl}\delta^{lm}+\hat{\omega}^j \hat{\omega}^l \hat{\omega}^l \hat{\omega}^m-\delta^{jl}\hat{\omega}^l \hat{\omega}^m-\delta^{lm}\hat{\omega}^j \hat{\omega}^l \\ =\delta^{jm}-\hat{\omega}^j \hat{\omega}^m$$
Best Answer
You don't need any series expansions to do this. $\mathbf Y^{-1}=\mathbf R(t+\Delta t) \mathbf R^{-1}(t)$ is a rotation matrix that rotates the body from the position at time $t$ to the position at time $t+\Delta t$. This is a rotation around the axis along $\boldsymbol\omega$ through the angle $|\boldsymbol\omega|\Delta t$. The trace of a rotation matrix with angle $\phi$ is $1+2\cos\phi$, so you can calculate $|\boldsymbol\omega|$ directly from the trace. To get the direction, you can solve the homogoeneous linear system $\mathbf Yx=x$.