I have 2 axis systems. One is just the natural axis system:
$e_{x} = \{1,0,0\}$, $e_{y} = \{0,1,0\}$, $e_{z} = \{0,0,1\}$
The other axis system is defined by 3 orthogonal vectors: $a$, $b$ and $c$.
What is the simplest method to find the Euler angles for this conversion, such that I can convert between axis systems with using a rotation matrix?
You don't need the Euler angles to obtain the rotation matrix (and its inverse) between the two systems. (There is no reason to work with Euler angles; they will only confuse you. Versors, or unit quaternions, are both better and easier to use to describe orientations.)
Let's say the other coordinate system basis — the three orthogonal unit vectors corresponding to the axes — are $$\begin{aligned} \hat{u} &= ( x_u , y_u , z_u ) \\ \hat{v} &= ( x_v , y_v , z_v ) \\ \hat{w} &= ( x_w , y_w , z_w ) \end{aligned}$$ Note that these are unit vectors, i.e. $\left\lVert \hat{u} \right\rVert = 1$, $\left\lVert\hat{v}\right\rVert = 1$, $\left\lVert\hat{w}\right\rVert = 1$. If your orthogonal vectors are of arbitrary (but nonzero) length, then simply normalize them. If you have some vector $\vec{a}$, $$\vec{a} = \left [ \begin{matrix} x_a \\ y_a \\ z_a \end{matrix} \right ] = ( x_a , y_a , z_a )$$ divide its components by its length, to normalize them: $$\hat{a} = \frac{\vec{a}}{\left\lVert\vec{a}\right\rVert} = \frac{\vec{a}}{\sqrt{\vec{a}\cdot\vec{a}}} = \left [ \begin{matrix} \frac{x_a}{\sqrt{x_a^2 + y_a^2 + z_a^2}} \\ \frac{y_a}{\sqrt{x_a^2 + y_a^2 + z_a^2}} \\ \frac{z_a}{\sqrt{x_a^2 + y_a^2 + z_a^2}} \end{matrix} \right ]$$
The rotation matrix $\mathbf{R}$ to rotate $\vec{p}$ in the natural coordinate system to $\vec{q}$ in the other coordinate system, via $\vec{q} = \mathbf{R} \vec{p}$, is $$\mathbf{R} = \left [ \begin{matrix} \hat{u} & \hat{v} & \hat{w} \end{matrix} \right ] = \left [ \begin{matrix} x_u & x_v & x_w \\ y_u & y_v & y_w \\ z_u & z_v & z_w \end{matrix} \right ]$$
As a pure rotation matrix, $\mathbf{R}$ is an orthonormal matrix, and its inverse is its transpose, $$\mathbf{R}^{-1} = \mathbf{R}^T = \left [ \begin{matrix} x_u & y_u & z_u \\ x_v & y_v & z_v \\ x_w & y_w & z_w \end{matrix} \right ]$$ This latter matrix rotates $\vec{q}$ back to $\vec{p}$ in the original coordinate system. In other words: $$\vec{q} = \mathbf{R} \vec{p} \quad \iff \quad \vec{p} = \mathbf{R}^T \vec{q}$$