ODE 45 / system of equations

Hi, I have just recently started using matlab and I am attempting to use ODE 45 to solve a system of non-linear ODEs.

I am having trouble with solving a system of equations so that I can run my function.

a truncated version of my main file looks like this:

clc; close all; clear all;

% Non-linear attitude equations
Tx=hd_vec(1)+hwxd+(obo_vec(2)*h_vec(3)-obo_vec(3)*h_vec(2))+(obo_vec(2)*hwz-obo_vec(3))*hwy==0;
Ty=hd_vec(2)+hwyd+(obo_vec(3)*h_vec(1)-obo_vec(1)*h_vec(3))+(obo_vec(3)*hwx-obo_vec(1))*hwz==0;
Tz=hd_vec(3)+hwzd+(obo_vec(1)*h_vec(2)-obo_vec(2)*h_vec(1))+(obo_vec(1)*hwy-obo_vec(2))*hwx==0;

% Variable isolation
phidd_Tx = isolate(Tx,phidd);
thetadd_Ty= isolate(Ty,thetadd);
psidd_Tz= isolate(Tz,psidd);

% Vectorization of inputs
x = [phi;phid;theta;thetad;psi;psid];

timerange=[0 120];%minutes
initialvalues=[1 0 1 0 1 0];
[t,x]=ode45(@(t,x) odefun_B_Sim1(t,x),timerange,initialvalues);
plot(t,y(:,1))
ylabel('test')
xlabel('test')

This code calls the function "odefun_B_Sim1" which is as follows:

function y=odefun_B_Sim1(t,x)

% Setup Variables
phi=x(1);
phid=x(2);
theta=x(3);
thetad=x(4);
psi=x(5);
psid=x(6);

% Setup Constants
Ixx=6; Iyy=8; Izz=4; Ixy=0; Ixz=0; Iyz=0; mu=(3.986*(10^14)); R=500;

% Setup Outputs
y1 = phid
y2 = -((Izz*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) - Iyz*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) + Ixz*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)))*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) - (Iyy*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) - Iyz*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) + Ixy*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)))*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) - Ixx*(sin(theta)*psidd + cos(theta)*psid*thetad + cos(phi)*cos(theta)*(mu/R^3)^(1/2)*phid - sin(phi)*sin(theta)*(mu/R^3)^(1/2)*thetad) + Ixz*((mu/R^3)^(1/2)*(sin(phi)*sin(psi)*psid - cos(phi)*cos(psi)*phid + cos(phi)*cos(psi)*sin(theta)*psid + cos(phi)*cos(theta)*sin(psi)*thetad - sin(phi)*sin(psi)*sin(theta)*phid) + sin(phi)*thetadd - cos(phi)*cos(theta)*psidd + cos(phi)*phid*thetad + cos(theta)*sin(phi)*phid*psid + cos(phi)*sin(theta)*psid*thetad) + hwy*(hwz*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) + sin(phi)*thetad - (mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - cos(phi)*cos(theta)*psid) + Ixy*((mu/R^3)^(1/2)*(cos(phi)*sin(psi)*sin(theta)*phid - cos(phi)*sin(psi)*psid - cos(psi)*sin(phi)*phid + cos(psi)*sin(phi)*sin(theta)*psid + cos(theta)*sin(phi)*sin(psi)*thetad) - cos(phi)*thetadd - cos(theta)*sin(phi)*psidd + sin(phi)*phid*thetad - cos(phi)*cos(theta)*phid*psid + sin(phi)*sin(theta)*psid*thetad) + 1)/Ixx
y3 = thetad
y4 = -((Izz*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) - Iyz*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) + Ixz*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)))*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)) - (Ixy*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) + Ixz*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) + Ixx*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)))*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) + Iyz*((mu/R^3)^(1/2)*(sin(phi)*sin(psi)*psid - cos(phi)*cos(psi)*phid + cos(phi)*cos(psi)*sin(theta)*psid + cos(phi)*cos(theta)*sin(psi)*thetad - sin(phi)*sin(psi)*sin(theta)*phid) - cos(phi)*cos(theta)*psidd + cos(phi)*phid*thetad + cos(theta)*sin(phi)*phid*psid + cos(phi)*sin(theta)*psid*thetad) - Iyy*((mu/R^3)^(1/2)*(cos(phi)*sin(psi)*sin(theta)*phid - cos(phi)*sin(psi)*psid - cos(psi)*sin(phi)*phid + cos(psi)*sin(phi)*sin(theta)*psid + cos(theta)*sin(phi)*sin(psi)*thetad) - cos(theta)*sin(phi)*psidd + sin(phi)*phid*thetad - cos(phi)*cos(theta)*phid*psid + sin(phi)*sin(theta)*psid*thetad) + hwz*(hwx*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) + sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)) + Ixy*(sin(theta)*psidd + cos(theta)*psid*thetad - phidd + cos(phi)*cos(theta)*(mu/R^3)^(1/2)*phid - sin(phi)*sin(theta)*(mu/R^3)^(1/2)*thetad) + 1)/(Iyz*sin(phi) + Iyy*cos(phi))
y5 = psid
y6 = -(Ixz*(cos(theta)*psid*thetad - phidd + cos(phi)*cos(theta)*(mu/R^3)^(1/2)*phid - sin(phi)*sin(theta)*(mu/R^3)^(1/2)*thetad) + (Ixy*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) + Ixz*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) + Ixx*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)))*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) - (Iyy*(cos(phi)*thetad - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) - Iyz*((mu/R^3)^(1/2)*(cos(psi)*sin(phi) - cos(phi)*sin(psi)*sin(theta)) - sin(phi)*thetad + cos(phi)*cos(theta)*psid) + Ixy*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)))*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)) - Izz*((mu/R^3)^(1/2)*(sin(phi)*sin(psi)*psid - cos(phi)*cos(psi)*phid + cos(phi)*cos(psi)*sin(theta)*psid + cos(phi)*cos(theta)*sin(psi)*thetad - sin(phi)*sin(psi)*sin(theta)*phid) + sin(phi)*thetadd + cos(phi)*phid*thetad + cos(theta)*sin(phi)*phid*psid + cos(phi)*sin(theta)*psid*thetad) + Iyz*((mu/R^3)^(1/2)*(cos(phi)*sin(psi)*sin(theta)*phid - cos(phi)*sin(psi)*psid - cos(psi)*sin(phi)*phid + cos(psi)*sin(phi)*sin(theta)*psid + cos(theta)*sin(phi)*sin(psi)*thetad) - cos(phi)*thetadd + sin(phi)*phid*thetad - cos(phi)*cos(theta)*phid*psid + sin(phi)*sin(theta)*psid*thetad) - hwx*(cos(phi)*thetad + hwy*(sin(theta)*psid - phid + cos(theta)*sin(phi)*(mu/R^3)^(1/2)) - (mu/R^3)^(1/2)*(cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(theta)) + cos(theta)*sin(phi)*psid) + 1)/(Ixz*sin(theta) + Izz*cos(phi)*cos(theta) - Iyz*cos(theta)*sin(phi))
y = [y1;y2;y3;y4;y5;y6]

The problem that I am getting is that y2 = phidd, y4 = thetadd, and y6 = psidd. This is problematic because y2 contains thetadd & psidd,

y4 contains phidd & psidd, and y6 contains phidd & thetadd. With 3 equations and 3 unknowns I should be able to solve it, but I cannot figure out how.

My files are attached.