How to solve the coupled differential equations with 3 variables? I have the following set of differential equations,
\begin{aligned}
& x' - \alpha(yz' - zy') =  -\beta(yP_z -zP_y) \\
& y' - \alpha(zx' - xz') =  -\beta(zP_x -xP_z) \\
& z' - \alpha(xy' - yx') =  -\beta(xP_y -yP_x)
\end{aligned}
I know how to implement runga kutta 4th order or Euler or Heun's method...what I can't figure out is how to update the $x'$ or $y'$ or $z'$ term while I am solving it. Any help or direction towards possible resources would be appreciated.
 A: Let $\vec{r} = (x,y,z)^T$ and $\hat{r}$ be corresponding unit vector.
Let $\vec{P} = (P_x,P_y,P_z)^T$ and $\vec{Q} = \vec{P} - (\vec{P}\cdot\hat{r})\hat{r}$. The ODE at hand can be rewritten as
$$\vec{r}' = \vec{r}\times( \alpha \vec{r}' - \beta \vec{P})
= \vec{r}\times(\alpha \vec{r}' - \beta \vec{Q})\tag{*1}
$$
Taking dot product with $\vec{r}$ on both sides, we find
$$\vec{r}\cdot\vec{r}' = 0 \implies |\vec{r}|^2 = R^2$$
for some constant $R$. Taking cross product with $\vec{r}$ on both sides of $(*1)$, we get
$$\require{cancel}
\begin{align}
\vec{r} \times \vec{r}' &= \vec{r}\times(\vec{r} \times(\alpha \vec{r}' - \beta \vec{Q}))\\
&= (\color{red}{\cancelto{0}{\color{gray}{\vec{r}\cdot(\alpha \vec{r}' - \beta \vec{Q})}}})\vec{r} - R^2 (\alpha \vec{r}' - \beta \vec{Q})\\
&= -R^2(\alpha \vec{r}' - \beta \vec{Q})
\end{align}
$$
Substitute this back into RHS of ($*1$) , we get
$$\begin{align}
\vec{r}' &= -\alpha R^2(\alpha \vec{r}' - \beta \vec{Q}) - \beta \vec{r} \times \vec{Q}\\
\iff \vec{r}' &=
\frac{\beta}{1 + \alpha^2R^2}(\alpha R^2 \vec{Q} - \vec{r} \times \vec{Q})
\end{align}
$$
In this form, the definition of $\vec{r}'$ is free of itself and you can throw it to any ODE solver you like.
A: On the left side you have a linear system in the first derivatives,
$$
\begin{bmatrix}
1&αz&-αy\\
-αz&1&αx\\
αy&-αx&1\\
\end{bmatrix}
\begin{bmatrix}
x'\\y'\\z'
\end{bmatrix}
=
\begin{bmatrix}
 -β(yP_z -zP_y)\\ -β(zP_x -xP_z) \\-β(xP_y -yP_x)
\end{bmatrix}
$$
In the evaluation of the ODE function, that is, the right side of the explicit first order system, you just have to construct these matrix and vector and call a solver for linear systems.
Then proceed as usual for explicit first-order systems.
