How to graph coupled differential equations

I have the following set of equations of motions for 2 masses and I am trying to plot the x displacement for each of the masses.

The MATLAB code i have for it is:

% Parameters
m1 = 4.0;
m2 = 6.0;
L = 1.5;
k = 100.0;  % Added missing parameter
g = 9.81;
F0 = 100;
tF = 1;

% Applied force F(t)
F = @(t) (t <= tF) * F0 .* sin(2*pi*t / tF);

% ODE function
odefun = @(t, y) [
    y(2);  % dx/dt
    compute_ddx(t, y, m1, m2, k, L, g, F);  % d²x/dt²
    y(4);  % dtheta/dt
    compute_ddtheta(t, y, m1, m2, k, L, g, F);  % d²theta/dt²
    y(6);  % dr/dt
    compute_ddr(t, y, m1, m2, k, L, g, F)  % d²r/dt²
    ];

% Initial state: [x, dx, theta, dtheta, r, dr]
y0 = [0; 0; 0; 0; L; 0];

% Duration
tspan = [0 3]; 

% Solve with tighter tolerances
options = odeset('RelTol',1e-6,'AbsTol',1e-8);
[t, Y] = ode45(odefun, tspan, y0, options);

% Extract states
x = Y(:,1); % x-displacement of m1 
theta = Y(:,3);
r = Y(:,5); % Use actual r(t) instead of constant L

% x displacement of m2 (x_m1 + r*sinθ)
x_m2 = x + r .* sin(theta);

% Create two separate figures
figure(1);
plot(x, t, 'b', 'LineWidth', 2)
xlabel('x m1 [m]');
ylabel('t [s]');
title('m1 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse');  
grid on;

figure(2);
plot(x_m2, t, 'r', 'LineWidth', 2)
xlabel('x m2 [m]');
ylabel('t [s]');
title('m2 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse'); 
grid on;

% Second derivative functions
function ddx = compute_ddx(t, y, m1, m2, k, L, g, F)
    Ft = F(t);
    theta = y(3); dtheta = y(4);
    r = y(5); dr = y(6);
    sinT = sin(theta); cosT = cos(theta);
    
    % % Corrected numerator with proper operators
    % numx = Ft - m2*(r*dtheta^2 + g*cosT + k*(L-r))*sinT - 2*m2*dr*dtheta*cosT + ...
    %        m2*r*(dtheta^2*sinT - (-2*dr*dtheta - g*sinT)*cosT/r);
    % 
    % % Corrected denominator
    % denx = m1 + m2 - m2*sinT^2 - m2*cosT^2;
    % ddx = numx/denx;
    ddx = (Ft - m2*k*(L-r)*sinT - 2*m2*r*dtheta^2*sinT)/( m1 + m2);
end

function ddtheta = compute_ddtheta(t, y, m1, m2, k, L, g, F)
    theta = y(3);
    r = y(5); dr = y(6);
    ddx = compute_ddx(t, y, m1, m2, k, L, g, F);
    dtheta = y(4);
    
    % Corrected ddtheta calculation
    ddtheta = (-ddx*cos(theta) - 2*dr*dtheta - g*sin(theta));
    ddtheta = ddtheta/r;
end

function ddr = compute_ddr(t, y, m1, m2, k, L, g, F)
    theta = y(3); dtheta = y(4);
    r = y(5);
    ddx = compute_ddx(t, y, m1, m2, k, L, g, F);
    
    % Corrected ddr calculation
    ddr = -ddx*sin(theta) + r*dtheta^2 + g*cos(theta) - k*(L-r)/m2;
end

Expectation

The graphs are supposed to come out as the green curves on these two graphs yet my code does not do that. I have checked that ddx, ddr and ddtheta have been computed correctly multiple times so the only thing i can think of it the way i am coding it.

Can someone pls assist?

Answer 1

[Torsten]

% Parameters
m1 = 4.0;
m2 = 6.0;
L = 1.5;
k = 100.0;  % Added missing parameter
g = 9.81;
F0 = 100;
tF = 1;
% Initial state: [x, dx, theta, dtheta, r, dr]
y0 = [0; 0; 0; 0; L; 0 ];
% Duration
tspan = [0 3]; 
% Solve with tighter tolerances
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t, Y] = ode45(@(t,y)odefun(t,y,m1, m2, k, L, g, F0, tF), tspan, y0, options);
% Extract states
x = Y(:,1); 
theta = Y(:,3);
r = Y(:,5); 
x_m2 = x + r .* sin(theta);
% Create two separate figures
figure(1);
plot(x, t, 'b', 'LineWidth', 2)
xlabel('x m1 [m]');
ylabel('t [s]');
title('m1 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse');  
grid on;

figure(2);
plot(x_m2, t, 'r', 'LineWidth', 2)
xlabel('x m2 [m]');
ylabel('t [s]');
title('m2 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse'); 
grid on

function dydt = odefun(t,y, m1, m2, k, L, g, F0, tF)
  x = y(1);
  xdot = y(2);
  theta = y(3);
  thetadot = y(4);
  r = y(5);
  rdot = y(6);

  F = (t <= tF) * F0 .* sin(2*pi*t / tF);
  
  A = zeros(3);
  b = zeros(3,1);
  A(1,1) = m1+m2;
  A(1,2) = m2*r*cos(theta);
  A(1,3) = m2*sin(theta);
  A(2,1) = m2*r*cos(theta);
  A(2,2) = m2*r^2;
  A(2,3) = 0;
  A(3,1) = m2*sin(theta);
  A(3,2) = 0;
  A(3,3) = m2;
  b(1) = F - 2*m2*rdot*thetadot*cos(theta) + m2*r*thetadot^2*sin(theta);
  b(2) = -2*m2*r*rdot*thetadot - m2*g*r*sin(theta);
  b(3) = m2*r*thetadot^2 + m2*g*cos(theta) + m2*k*(L-r);
  ddot = A\b;
  
  xddot = ddot(1);
  thetaddot = ddot(2);
  rddot = ddot(3);
  
  dydt = [xdot;xddot;thetadot;thetaddot;rdot;rddot];
end