MATLAB: Using ode45 to model equations of motion for a coupled nonlinear system of ordinary differential equations

ode45

I am trying to model these equations of motions (EoMs) shown in the attached picture using MATLAB's ode45 function. The responses that I am interested in modeling are and , all other variables are constants in the equations. Since it is a second order differential equation, I convert the system of equations from 2nd order to 1st order in order to model the EoMs. However, when I run my code I do not get a response in the direction like I am anticipating that I should since the EoMs are coupled. This is shown in the 2nd figure attached.
I set and solve for , by plugging in the EoMs and using MATLAB's solve function to get , in terms of . I've also attached my MATLAB code. Any help would be greatly appreciated. I've been stuck on this for a while and it's starting to drive me bonkers so a second set of eyes or opinion on how to approach the problem would help. 
clear all; close all 
% Solve for xdot2 and xdot4 in terms of x1, x2, x3, and x4
syms xdot2 omegas zetas x1 x2 x3 x4 xdot4 cs OMEG X0 R t alphas n
eqn1 = xdot2 + 2*omegas*zetas*(x2 + cs^2*(x1^2*x2 + x1*x3*x4)) + cs^2*(x1^2*xdot2 + x1*x2^2 + x1*x3*xdot4 + x1*x4^2) + ...
    omegas^2*(1 + alphas*(x1^2 + x3^2)^(n-1))*x1 - OMEG^2*(X0/R)*cos(OMEG*t) == 0; 
eqn2 = xdot4 + 2*omegas*zetas*(x4 + cs^2*(x3^2*x4 + x1*x2*x3)) + cs^2*(x3^2*xdot4 + x3*x4^2 + x1*xdot2*x3 + x2^2*x3) + ...
    omegas^2*(1 + alphas*(x1^2 + x3^2)^(n-1))*x3 == 0; 
xdots = solve([eqn1, eqn2],[xdot4, xdot2]);
% Plug xdot2 and xdot4 into the bauer function defined below
xdot2 = simplify(xdots.xdot2); 
xdot4 = simplify(xdots.xdot4); 
 
% initial conditions 
tspan = 0:.025:2*pi; 
z0 = [0 0 0 0.0]; 
[t,z] = ode45(@bauer, tspan, z0); 
% plot results 
figure()
subplot(211)
plot(t,z(:,1),'+'); 
grid on 
xlabel('Time (sec)'); ylabel('\eta_s'); 
subplot(212)
plot(t,z(:,3)); 
grid on 
xlabel('Time (sec)'); ylabel('\zeta_s'); 
function EoM = bauer(t,x)
omegas = 10; 
zetas = .1; 
OMEG = 11; 
n = 2; 
alphas =  2/3; 
R = 5.938; 
X0 = 0.032;
cs = 2; 
x1 = x(1); 
x2 = x(2); 
x3 = x(3); 
x4 = x(4); 
xd2 = -(R*omegas^2*x1 - OMEG^2*X0*cos(OMEG*t) + 2*R*omegas*x2*zetas + R*cs^2*x1*x2^2 + R*cs^2*x1*x4^2 - OMEG^2*X0*cs^2*x3^2*cos(OMEG*t) + R*alphas*omegas^2*x1*(x1^2 + x3^2)^(n - 1) + 2*R*cs^2*omegas*x1^2*x2*zetas + 2*R*cs^2*omegas*x2*x3^2*zetas)/(R*(cs^2*x1^2 + cs^2*x3^2 + 1));
xd4 = -(R*omegas^2*x3 + 2*R*omegas*x4*zetas + R*cs^2*x2^2*x3 + R*cs^2*x3*x4^2 + R*alphas*omegas^2*x3*(x1^2 + x3^2)^(n - 1) + OMEG^2*X0*cs^2*x1*x3*cos(OMEG*t) + 2*R*cs^2*omegas*x1^2*x4*zetas + 2*R*cs^2*omegas*x3^2*x4*zetas)/(R*(cs^2*x1^2 + cs^2*x3^2 + 1)); 
 
EoM = [x2; xd2; x4; xd4]; 
end

Best Answer

subplot(212)
R = 5.938;
plot(t,z(:,3)*R); 
grid on 
xlabel('Time (sec)'); ylabel('\zeta_s');
Related Question