MATLAB: Ode45 problems — not getting expected solution

Hello,

I have a set of nonlinear satellite attitude equations to solve. I've implemented this as I've done in the past (at least I think I've done everything the same), and the solutions I get are linear. There should be an oscillation in each of the angles phi, tht, psi over time. My suspicion is that I didn't define the vector correctly being sent to ode45, my ode function doesn't update the contents of this vector correctly or there's some sort of scaling issue. I've asked other classmates already, the TA for my class, and stared at it for hours without seeing the error, so I'm turning to the experts for assistance.

Thank you,

Rick

% Derive Equations of Motion
clc; clear;
% PRY = [phi,tht,psi,dphi,dtht,dpsi]';
% Create sym variables
syms Sphi(t) Stht(t) Spsi(t) dSphi(t) dStht(t) dSpsi(t) ddSphi(t) ddStht(t)...
    ddSpsi(t)
assume([Sphi(t) Stht(t) Spsi(t) dSphi(t) dStht(t) dSpsi(t) ddSphi(t) ...
    ddStht(t) ddSpsi(t)],'real')
syms SOmega SR SG SmE SIxx SIyy SIzz real positive
% Form rotation matrices and cosine direction matrices
Rxphi = [1,0,0;0,cos(Sphi),sin(Sphi);0,-sin(Sphi),cos(Sphi)];
Rytht(t) = [cos(Stht),0,-sin(Stht);0,1,0;sin(Stht),0,cos(Stht)];
Rzpsi(t) = [cos(Spsi),sin(Spsi),0;-sin(Spsi),cos(Spsi),0;0,0,1];
BCA = Rxphi*Rytht*Rzpsi;
ACB = transpose(BCA);
% Compute O_omega_A in B-frame
OwA_B = BCA*[0;0;SOmega];
% Compute A_omega_B in B-frame
AwBcross = BCA*diff(ACB,t);
AwB_Bx = [0,0,1]*AwBcross*[0;1;0];
AwB_By = [1,0,0]*AwBcross*[0;0;1];
AwB_Bz = [0,1,0]*AwBcross*[1;0;0];
AwB_B = [ AwB_Bx;AwB_By;AwB_Bz ];
% Compute O_omega_B in B-frame
OwB_B = OwA_B + AwB_B;
% Decompose O_omega_B
OwB = eye(3)*OwB_B(t);
OwB = [OwB(1,:);OwB(2,:);OwB(3,:)];
% Moment of Inertia Matrix
I = [SIxx,0,0;0,SIyy,0;0,0,SIzz];
% Angular Momentum
h = I*OwB;
hdot = diff(h,t);
% Torques
K = 3*SG*SmE/SR^5;
Gx = K*(sin(2*Sphi)*(cos(Stht))^2*(SIzz-SIyy));
Gy = K*(sin(2*Stht)*cos(Sphi)*(SIzz-SIxx));
Gz = K*(sin(2*Stht)*sin(Sphi)*(SIxx-SIyy));
Tau = hdot + cross(OwB,h);
TauX = Tau(1,:);
TauY = Tau(2,:);
TauZ = Tau(3,:);
% Substitute
TauX = subs(TauX,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauY = subs(TauY,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauZ = subs(TauZ,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauX = subs(TauX,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauY = subs(TauY,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauZ = subs(TauZ,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauX = subs(TauX,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
TauY = subs(TauY,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
TauZ = subs(TauZ,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
% Tx = TauX == Gx;
% Ty = TauY == Gy;
% Tz = TauZ == Gz;
Tx = TauX == 0;
Ty = TauY == 0;
Tz = TauZ == 0;
isolate(Tx,ddSphi)
isolate(Ty,ddStht)
isolate(Tz,ddSpsi)



% Main
% Sets intitial values, creates PRY vector of unknowns and paramters; calls
% other functions and solver; sorts data and plots results

clc; clear; close all

% Set initial values

PRY = zeros(6,1);

phi = 0; tht = 0; psi = 0;
dphi = 0; dtht = 0; dpsi = 0;
RE = 6378.137e3 ;                       % m - radius of earth
RS = 500e3;                             % m - orbital altitude
R = (RE+RS);                            % m
mE = 5.972e24;                          % kg
G = 6.6743e-11;                         % m^3/kg s^2
Omega = sqrt(G*mE/R^3);                 % rad/s
Ixx = 6; Iyy = 8; Izz = 4;              % kg m^2
mu = G*mE;
per = 2*pi*(sqrt(R^3/mu));
tend = 1.6*per;
tspan = [0 tend];
ic = [1;5;10;15;20;25];

for k= 1:6
    phi =ic(k);
    tht=ic(k);
    psi=ic(k);
    PRY = [phi,tht,psi,dphi,dtht,dpsi]';
    [t,prydot] = ode45(@(t,pry) proj_ODE(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
%     [t,prydot] = ode45(@(t,pry) proj_ODE_L(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
    
    % Plotting results; pitch, roll, yaw for each initial condition
    xticks = linspace(0,1.6,9);
    
    figure(k)
    %     prydot = rad2deg(prydot);
    subplot(3,1,1)
    plot(t,prydot(:,1))
    grid on
    title('Pitch, \phi');
    xlabel('Time, [s]'); ylabel('');
    
    subplot(3,1,2)
    plot(t,prydot(:,2));
    grid on
    title('Roll, \theta');
    xlabel('Time, [s]');
    
    subplot(3,1,3)
    plot(t,prydot(:,3));
    grid on
    title('Yaw, \psi');
    xlabel('Time, [s]'); ylabel('Attitude, [rad]');
    
end



% pry_ODE.m
% myODE function for solver
% PRY = [phi,tht,psi,dphi,dtht,dpsi];

function prydot = proj_ODE(t,pry,Omega,R,mE,G,Ixx,Iyy,Izz)
    
    prydot = zeros(size(pry));
    
    phi = pry(1);
    tht = pry(2);
    psi = pry(3);
    dphi = pry(4);
    dtht = pry(5);
    dpsi = pry(6);
    
    K = 3*G*mE/R^5;
    w0 = Omega;
    
    Gx = K*(sin(2*phi)*(cos(tht))^2*(Izz-Iyy));
    Gy = K*(sin(2*tht)*cos(phi)*(Izz-Ixx));
    Gz = K*(sin(2*tht)*sin(phi)*(Ixx-Iyy));
    
    prydot(1) = dphi;
    prydot(2) = dtht;
    prydot(3) = dpsi;

    prydot(4) = (Ixx*((sin(phi)*sin(psi) ...
        + cos(phi)*cos(psi)*sin(tht))*(cos(phi)*cos(psi)*prydot(6) ...
        - cos(phi)*sin(psi)*dphi^2 - cos(phi)*sin(psi)*dpsi^2 ...
        + cos(psi)*sin(phi)*sin(tht)*dphi^2 + cos(psi)*sin(phi)*sin(tht)*dpsi^2 ...
        + cos(psi)*sin(phi)*sin(tht)*dtht^2 - 2*cos(psi)*sin(phi)*dphi*dpsi ...
        - cos(psi)*cos(tht)*sin(phi)*prydot(5) + sin(phi)*sin(psi)*sin(tht)*prydot(6) ...
        - 2*cos(phi)*cos(psi)*cos(tht)*dphi*dtht ...
        + 2*cos(phi)*sin(psi)*sin(tht)*dphi*dpsi ...
        + 2*cos(tht)*sin(phi)*sin(psi)*dpsi*dtht) - (cos(psi)*sin(phi) ...
        - cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*prydot(6) ...
        + cos(phi)*cos(psi)*dphi^2 + cos(phi)*cos(psi)*dpsi^2 ...
        + sin(phi)*sin(psi)*sin(tht)*dphi^2 + sin(phi)*sin(psi)*sin(tht)*dpsi^2 ...
        + sin(phi)*sin(psi)*sin(tht)*dtht^2 - 2*sin(phi)*sin(psi)*dphi*dpsi ...
        - cos(psi)*sin(phi)*sin(tht)*prydot(6) - cos(tht)*sin(phi)*sin(psi)*prydot(5) ...
        - 2*cos(phi)*cos(psi)*sin(tht)*dphi*dpsi - 2*cos(phi)*cos(tht)*sin(psi)*dphi*dtht ...
        - 2*cos(psi)*cos(tht)*sin(phi)*dpsi*dtht) - (cos(phi)*sin(psi)*dphi ...
        + cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
        - cos(psi)*sin(phi)*sin(tht)*dphi ...
        - cos(phi)*sin(psi)*sin(tht)*dpsi)*(sin(phi)*sin(psi)*dphi ...
        - cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
        + cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
        - (sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
        + cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
        - sin(phi)*sin(psi)*sin(tht)*dphi)*(cos(phi)*sin(psi)*sin(tht)*dphi ...
        - cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
        + cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
        + cos(phi)*cos(tht)*(sin(phi)*sin(tht)*prydot(5) ...
        + cos(tht)*sin(phi)*dphi^2 + cos(tht)*sin(phi)*dtht^2 ...
        + 2*cos(phi)*sin(tht)*dphi*dtht) + w0*cos(tht)*dtht ...
        + cos(tht)*sin(phi)*dphi*(cos(phi)*cos(tht)*dphi - sin(phi)*sin(tht)*dtht) ...
        + cos(phi)*sin(tht)*dtht*(cos(phi)*cos(tht)*dphi - sin(phi)*sin(tht)*dtht)) ...
        + Iyy*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
        + w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
        + cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
        - cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
        + cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
        + cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
        - sin(phi)*sin(psi)*sin(tht)*dphi))*((cos(phi)*cos(psi) ...
        + sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
        - sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
        - cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
        + cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) - cos(tht)^2*sin(phi)*dtht) ...
        - Izz*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
        + w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
        + cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
        - cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
        + cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
        + cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
        - sin(phi)*sin(psi)*sin(tht)*dphi))*((cos(phi)*cos(psi) ...
        + sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
        - sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
        - cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
        + cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
        - cos(tht)^2*sin(phi)*dtht))/(Ixx*(cos(phi)^2*cos(tht)^2 + (sin(phi)*sin(psi) ...
        + cos(phi)*cos(psi)*sin(tht))^2 + (cos(psi)*sin(phi) ...
        - cos(phi)*sin(psi)*sin(tht))^2));
    
  prydot(5) = (Iyy*(cos(tht)*sin(psi)*(cos(phi)*cos(psi)*prydot(4) ...
        - sin(phi)*sin(psi)*prydot(6) - cos(psi)*sin(phi)*dphi^2 ...
        - cos(psi)*sin(phi)*dpsi^2 + cos(phi)*sin(psi)*sin(tht)*dphi^2 ...
        + cos(phi)*sin(psi)*sin(tht)*dpsi^2 + cos(phi)*sin(psi)*sin(tht)*dtht^2 ...
        - 2*cos(phi)*sin(psi)*dphi*dpsi - cos(phi)*cos(psi)*sin(tht)*prydot(6) ...
        + sin(phi)*sin(psi)*sin(tht)*prydot(4) ...
        - 2*cos(phi)*cos(psi)*cos(tht)*dpsi*dtht ...
        + 2*cos(psi)*sin(phi)*sin(tht)*dphi*dpsi ...
        + 2*cos(tht)*sin(phi)*sin(psi)*dphi*dtht) ...
        - sin(tht)*(cos(tht)*sin(phi)*prydot(4) + cos(phi)*cos(tht)*dphi^2 ...
        + cos(phi)*cos(tht)*dtht^2 - 2*sin(phi)*sin(tht)*dphi*dtht) ...
        + cos(psi)*cos(tht)*(sin(phi)*sin(psi)*dphi^2 + sin(phi)*sin(psi)*dpsi^2 ...
        - cos(phi)*sin(psi)*prydot(4) - cos(psi)*sin(phi)*prydot(6) ...
        + cos(phi)*cos(psi)*sin(tht)*dphi^2 + cos(phi)*cos(psi)*sin(tht)*dpsi^2 ...
        + cos(phi)*cos(psi)*sin(tht)*dtht^2 - 2*cos(phi)*cos(psi)*dphi*dpsi ...
        + cos(psi)*sin(phi)*sin(tht)*prydot(4) ...
        + cos(phi)*sin(psi)*sin(tht)*prydot(6) ...
        + 2*cos(psi)*cos(tht)*sin(phi)*dphi*dtht ...
        + 2*cos(phi)*cos(tht)*sin(psi)*dpsi*dtht ...
        - 2*sin(phi)*sin(psi)*sin(tht)*dphi*dpsi) ...
        - cos(tht)*dtht*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
        + w0*sin(phi)*sin(tht)*dtht - cos(psi)*cos(tht)*dpsi*(sin(phi)*sin(psi)*dpsi ...
        - cos(phi)*cos(psi)*dphi + cos(phi)*cos(psi)*sin(tht)*dpsi ...
        + cos(phi)*cos(tht)*sin(psi)*dtht - sin(phi)*sin(psi)*sin(tht)*dphi) ...
        + cos(tht)*sin(psi)*dpsi*(cos(phi)*sin(psi)*dphi + cos(psi)*sin(phi)*dpsi ...
        + cos(phi)*cos(psi)*cos(tht)*dtht - cos(psi)*sin(phi)*sin(tht)*dphi ...
        - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
        + cos(psi)*sin(tht)*dtht*(cos(phi)*sin(psi)*dphi ...
        + cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
        - cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
        + sin(psi)*sin(tht)*dtht*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
        + cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
        - sin(phi)*sin(psi)*sin(tht)*dphi) - w0*cos(phi)*cos(tht)*dphi) ...
        - Ixx*((cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
        - sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
        - cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
        + cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
        - cos(tht)^2*sin(phi)*dtht)*((sin(phi)*sin(psi) ...
        + cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
        - cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
        + cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
        - (cos(psi)*sin(phi) ...
        - cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
        - cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
        + cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
        - w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
        - sin(phi)*sin(tht)*dtht)) + Izz*((cos(phi)*cos(psi) ...
        + sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
        - sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
        - cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
        + cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
        - cos(tht)^2*sin(phi)*dtht)*((sin(phi)*sin(psi) ...
        + cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
        - cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
        + cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
        - (cos(psi)*sin(phi) ...
        - cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
        - cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
        + cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
        - w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
        - sin(phi)*sin(tht)*dtht)))/(Iyy*(cos(phi)*cos(psi)^2*cos(tht)^2 ...
        + cos(phi)*cos(tht)^2*sin(psi)^2 + cos(phi)*sin(tht)^2));
    
    prydot(6) = (Izz*((cos(phi)*cos(psi) ...
        + sin(phi)*sin(psi)*sin(tht))*(sin(psi)*sin(tht)*prydot(5) ...
        + cos(tht)*sin(psi)*dpsi^2 + cos(tht)*sin(psi)*dtht^2 ...
        + 2*cos(psi)*sin(tht)*dpsi*dtht) + (cos(tht)*sin(psi)*dpsi ...
        + cos(psi)*sin(tht)*dtht)*(sin(phi)*sin(psi)*dphi ...
        - cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
        + cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
        - (cos(psi)*cos(tht)*dpsi ...
        - sin(psi)*sin(tht)*dtht)*(cos(phi)*sin(psi)*sin(tht)*dphi ...
        - cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
        + cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
        - (cos(phi)*sin(psi) ...
        - cos(psi)*sin(phi)*sin(tht))*(cos(psi)*sin(tht)*prydot(5) ...
        + cos(psi)*cos(tht)*dpsi^2 + cos(psi)*cos(tht)*dtht^2 ...
        - 2*sin(psi)*sin(tht)*dpsi*dtht) + cos(tht)^2*sin(phi)*prydot(5) ...
        - 2*cos(tht)*sin(phi)*sin(tht)*dtht^2 + cos(phi)*cos(tht)^2*dphi*dtht ...
        + w0*cos(tht)*sin(phi)*dphi + w0*cos(phi)*sin(tht)*dtht) ...
        + Ixx*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
        + w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
        + cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
        - cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
        + cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
        + cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
        - sin(phi)*sin(psi)*sin(tht)*dphi))*((sin(phi)*sin(psi) ...
        + cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
        - cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
        + cos(psi)*cos(tht)*sin(phi)*d
            

            
            
            

    
        
	
    
        Best Answer
    

    
								
								    K = 3*G*mE/R^5;
It's been a while since I've studied any physics, but this looks odd to me.
RE = 6378.137e3 ;                       % m - radius of earth
RS = 500e3;                             % m - orbital altitude
R = (RE+RS);                            % m
mE = 5.972e24;                          % kg
G = 6.6743e-11;                         % m^3/kg s^2
By dimensional analysis, K has units of or  times whatever the unit of 3 is. If that's supposed to be the gravitational force between a 3 kg satellite and the Earth, shouldn't that be over R^2 (which would give Newtons, ) instead of over R^5?
						    
								
    
        
    

Related Solutions
MATLAB: Ode45 Error using odearguments (line 95)

								
								First, specify: 
prydot = zeros(6,1);
to create the necessary column vector.  
Then, since it is only neccessary for ‘pry_ODE’ to return the vector of derivatives, not the time as well, change the function declaration to: 
prydot = pry_ODE(t,pry,Omega,R,mE,G,Ixx,Iyy,Izz)
and your system runs without error.  
						    
						
MATLAB: Numerical Integration in matlab

								
								There is no reliable numeric integration method in any programming language. For any given numeric integration method, it is possible to construct a function which the numeric integration method will return an arbitrarily wrong solution. 
Trapazoidal Rule and Simpson's Rule are not reliable at all.
You could consider using integral(); provide AbsTol and RelTol parameters, and provide WayPoints whenever meaningful.
You could consider using vpaintegral() with similar parameters.
						    
						
        
        
            
            
            Related Question
Symbolic Derivative in matlab
Driven Damped Pendulum Axes Question
Find the derivative of a piecewise function and plot the result
Error:An array for multiple LHS assignment cannot contain expressions.
Solve and plot a system of nonlinear 2nd order differential equations