I have 2 coupled systems of differential equations, used for spacecraft simulation, which are to be integrated simultaneously:

In this system, omega_i and epsilon_j are the variables. (i belongs to {1,2,3} and j belongs to {1,2,3,4}). All other symbols are constants.

In MATLAB, I typically make a new vector as [w;e], which has both the w and e terms, and use it in a single function.

```
function dwdt = KDE_dwdt(t,w,I)
dwdt = zeros(7,1);
dwdt(1) = -w(3)*w(2)*(I(3)-I(2))/I(1);
dwdt(2) = -w(3)*w(1)*(I(1)-I(3))/I(2);
dwdt(3) = -w(1)*w(2)*(I(2)-I(1))/I(3);
e = w(4:7);
de = [e(4) -e(3) e(2) e(1);...
e(3) e(4) -e(1) e(2);...
-e(2) e(1) e(4) e(3);...
-e(1) -e(2) -e(3) e(4); ];
dwdt(4:7) = 1/2 * de * [w(1:3);0];
end
```

Should I follow a similar way in Julia. I don’t know if Julia allows for 2 separate functions to be integrated simultaneously, or if there are more efficient ways to solve this system. This is because I already have the individual system of equations written for omega, and epsilon, but I’m not sure if those can be solved with individually written functions.

For instance, in Julia, I have the following function written for epsilon and omegas:

```
function omega_propagator(dω,ω,t,I,L)
dω[1] = (-ω[3]*ω[2]*(I[3]-I[2]))/I[1];
dω[2] = (-ω[3]*ω[1]*(I[1]-I[3]))/I[2];
dω[3] = (-ω[1]*ω[2]*(I[2]-I[1]))/I[3];
end
function KDE_EP!(de,e,wf,t)
ω = wf(de,e,t); #here wf was a function which gave values of omega, but not integrated numerically
de.= 0.5*[e[4] -e[3] e[2] e[1];e[3] e[4] -e[1] e[2]; -e[2] e[1] e[4] e[3];-e[1] -e[2] -e[3] e[4]]*[ω;0];
end
```

Is there a way that I can integrate omega values in `KDE_EP!`

and use omega values in `omega_propagator`

for the next step?

Thanks!