We'have derived analytical equations to solve the motion for the system. Now we'll write code in Matlab by using Matlab Live Script. 

clear 

clc
close all 
warning('off')
syms t 
syms x_1(t) x_2(t) theta(t) X(t) Y(t) 
syms x_dot_1 x_dot_2 theta_dot X_dot Y_dot
x_dot_1   = diff(x_1, t);
x_dot_2   = diff(x_2, t);
theta_dot = diff(theta, t);
X_dot     = diff(X, t);
Y_dot     = diff(Y, t);
syms m_1 m_2 m_3 m_4 m_5 
syms I_3 I_4
syms k_1 k_2 k_3 c_1
syms l h r g 
syms mu N 
syms F(t) 
syms lambda_1 lambda_2

disp('--------------------KINETIC ENERGY--------------------')
K_1 = 1/2 * m_1 * x_dot_1^2
K_2 = 1/2 * m_2 * x_dot_1^2
K_3 = 1/2 * m_3 * x_dot_1^2  + 1/2 * I_3 * theta_dot^2
x_4 = x_1 + (r+l/2)*cos(theta)
y_4 = (h+l)+(r+l/2)*sin(theta)
v_4_square_ = diff(x_4,t)^2 + diff(y_4,t)^2
K_4 = (1/2 * m_4 * v_4_square_) + 1/2 * I_4 * theta_dot^2;
K_5 = 1/2 * m_5 * (X_dot^2 + Y_dot^2);
K = K_1 + K_2 + K_3 + K_4 + K_5 

disp('--------------------POTENTIAL ENERGY--------------------')
U_1 = m_1 * g * h/2 + (1/2 * k_1 * x_1^2)
U_2 = m_2 * g * (h + l/2) 
U_3 = m_3 * g * (h + l) + (1/2 * k_2 * theta^2)
U_4 = m_4 * g * (h + l + (r + l/2)*sin(theta))
U_5 = m_5 * g * Y + 1/2* k_3 * x_2^2 
U = U_1 + U_2 + U_3 + U_4 + U_5
disp('--------------------DISSIPATIVE ENERGY--------------------')
% No friction force
D = 1/2 * c_1 * x_dot_1^2 ; %+ sign(x_dot_2) * mu *  m_5*((r+x_2)*diff(theta_dot,t) + 2*x_dot_2*theta_dot + g*cos(theta))
disp('--------------------EXTERNAL FORCES--------------------')

Q_x_1 = 0   +  F * sin(theta)
Q_x_2 = 0  
Q_theta = F * (l+r) 
Q_X = 0 % - F * cos(theta) % ???
Q_Y = 0 % - F * sin(theta) % ???

disp('--------------------CONSTRAINTS --------------------')
phi_1 = X - (x_1 + (r+x_2)*cos(theta))
phi_2 = Y - (h + l + (r+x_2)*sin(theta)) 

d_phi_1_dx_1 = diff(phi_1,x_1)
d_phi_2_dx_1 = diff(phi_2,x_1)

d_phi_1_dx_2 = diff(phi_1,x_2)
d_phi_2_dx_2 = diff(phi_2,x_2)

d_phi_1_dtheta = diff(phi_1,theta)
d_phi_2_dtheta = diff(phi_2,theta)

d_phi_1_dX = diff(phi_1,X)
d_phi_2_dX = diff(phi_2,X)

d_phi_1_dY = diff(phi_1,Y)
d_phi_2_dY = diff(phi_2,Y)

PfaffianForm = [d_phi_1_dx_1, d_phi_1_dx_2, d_phi_1_dtheta, d_phi_1_dX, d_phi_1_dY; 
                d_phi_2_dx_1, d_phi_2_dx_2, d_phi_2_dtheta, d_phi_2_dX, d_phi_2_dY]

cons_1 = lambda_1 * d_phi_1_dx_1 + lambda_2 * d_phi_2_dx_1
cons_2 = lambda_1 * d_phi_1_dx_2 + lambda_2 * d_phi_2_dx_2
cons_3 = lambda_1 * d_phi_1_dtheta + lambda_2 * d_phi_2_dtheta
cons_4 = lambda_1 * d_phi_1_dX + lambda_2 * d_phi_2_dX
cons_5 = lambda_1 * d_phi_1_dY + lambda_2 * d_phi_2_dY




disp('--------------------DERIVE LAGRANGIAN--------------------')
% ---------------------
dK_d_x_dot_1 = diff(K,x_dot_1) % dK/dx_dot_1
d_dt_dK_d_x_dot_1 = diff(dK_d_x_dot_1,t) % d/dt(dK/dx_dot_1)
dK_d_x_1 = diff(K, x_1) % dK/d_x_1
dU_d_x_1 = diff(U, x_1) % dU/d_x_1
Dd_d_x_dot_1 = diff(D,x_dot_1) %dD/d_x_dot_1
L_1 = d_dt_dK_d_x_dot_1 - dK_d_x_1 + dU_d_x_1 + Dd_d_x_dot_1 == Q_x_1 + cons_1 
% ---------------------
dK_d_x_dot_2 = diff(K,x_dot_2); % dK/dx_dot_2
d_dt_dK_d_x_dot_2 = diff(dK_d_x_dot_2,t); % d/dt(dK/dx_dot_2)
dK_d_x_2 = diff(K, x_2); % dK/d_x_2
dU_d_x_2 = diff(U, x_2); % dU/d_x_2
Dd_d_x_dot_2 = diff(D,x_dot_2); %dD/d_x_dot_2
L_2 = d_dt_dK_d_x_dot_2 - dK_d_x_2 + dU_d_x_2 + Dd_d_x_dot_2 == Q_x_2 + cons_2
% ---------------------
dK_d_theta_dot = diff(K,theta_dot); % dK/dtheta_dot
d_dt_dK_d_theta_dot = diff(dK_d_theta_dot,t); % d/dt(dK/dtheta_dot)
dK_d_theta = diff(K, theta); % dK/d_theta
dU_d_theta = diff(U, theta); % dU/d_theta
Dd_d_theta_dot = diff(D,theta_dot); %dD/d_theta_dot
L_3 = d_dt_dK_d_theta_dot - dK_d_theta + dU_d_theta + Dd_d_theta_dot == Q_theta + cons_3
% ---------------------
dK_d_X_dot = diff(K,X_dot); % dK/dX_dot
d_dt_dK_d_X_dot = diff(dK_d_X_dot,t); % d/dt(dK/dX_dot)
dK_d_X = diff(K, X); % dK/d_X
dU_d_X = diff(U, X); % dU/d_X
Dd_d_X_dot = diff(D,X_dot); %dD/d_X_dot
L_4 = d_dt_dK_d_X_dot - dK_d_X + dU_d_X + Dd_d_X_dot == Q_X + cons_4
% ---------------------
dK_d_Y_dot = diff(K,Y_dot); % dK/dY_dot
d_dt_dK_d_Y_dot = diff(dK_d_Y_dot,t); % d/dt(dK/dY_dot)
dK_d_Y = diff(K, Y); % dK/d_Y
dU_d_Y = diff(U, Y); % dU/d_Y
Dd_d_Y_dot = diff(D,Y_dot); %dD/d_Y_dot
L_5 = d_dt_dK_d_Y_dot - dK_d_Y + dU_d_Y + Dd_d_Y_dot == Q_Y + cons_5 
% ---------------------
EQ = [L_1; L_2; L_3; L_4; L_5]
EQ_S = [EQ; phi_1 == 0; phi_2 == 0]
 
Find equilibrium 
% Find equilibrium 
EQ_N = EQ; 
EQ_N = subs(EQ_N, x_dot_1,0); 
EQ_N = subs(EQ_N, x_dot_2,0);
EQ_N = subs(EQ_N, theta_dot,0);
EQ_N = subs(EQ_N, X_dot,0);
EQ_N = subs(EQ_N, Y_dot,0);

EQ_N = subs(EQ_N, diff(x_dot_1,t),0); 
EQ_N = subs(EQ_N, diff(x_dot_2,t),0);
EQ_N = subs(EQ_N, diff(theta_dot,t),0);
EQ_N = subs(EQ_N, diff(X_dot,t),0);
EQ_N = subs(EQ_N, diff(Y_dot,t),0)


x_ddot_1 = diff(x_dot_1,t);
x_ddot_2 = diff(x_dot_2,t);
theta_ddot = diff(theta_dot,t);
vars = [x_1 x_dot_1 x_ddot_1 x_2 x_dot_2 x_ddot_2 theta theta_dot theta_ddot]; 
pars = [m_1 m_2 m_3 m_4 m_5 k_1 k_2 k_3 c_1 I_3 I_4 h l r g F mu]; 
CQ=[L_1; L_2; L_3];
X = (x_1 + (r+x_2)*cos(theta))
Y = (h + l + (r+x_2)*sin(theta)) 
lambda_1_new = m_5 * diff(X,t,2);
lambda_2_new = m_5 * (g + diff(Y,t,2));
CQ = subs(CQ,lambda_1,lambda_1_new);
CQ = subs(CQ,lambda_2,lambda_2_new)

 
CQ = subs(CQ, ...
    [m_1, m_2, m_3, m_4, m_5, k_1, k_2,        k_3,   c_1, I_3, I_4,        h,     l,   r,    g,], ...
    [8,   2,   6,   2,   3,   50,   250,       200,   10,   3,  8.6667,   0.5,     2,   1,    -9.81]);
[newEQS, newVARS] = reduceDifferentialOrder(CQ,[x_1(t) x_2(t) theta(t) F(t)]);
[M,zxc] = massMatrixForm(newEQS,newVARS);
f = M\zxc;
odefun = odeFunction(f,newVARS,'File','ode_solver.m');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_1_IC = 0; 
x_2_IC = 0; 
theta_IC = 0; 
F_IC = +100; 
x_dot_1_IC = 0; 
x_dot_2_IC = 0;
theta_dot_IC = 0 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initalConditions =   [x_1_IC,  x_2_IC, theta_IC, F_IC, x_dot_1_IC, x_dot_2_IC, theta_dot_IC]; 
tspan = [0 200]; 
[time_, result_] = ode45(odefun, tspan, initalConditions); 
result_x_1 = result_(:,1); 
result_x_2 = result_(:,2); 
result_theta = result_(:,3);
result_F = result_(:,4); 
result_x_dot_1 = result_(:,5); 
result_x_dot_2 = result_(:,6);
result_theta_dot = result_(:,7); 
clc
disp("-------------------------------------------------")
disp("-------------LAGRANGE FINISHED-------------------")
model_params
disp("-------------------------------------------------")
disp("-------------SIMULATION FINISHED-------------------")

close all 
figure(1)
subplot(3,1,1)
plot(out.simulation_results.Time(:), out.simulation_results.Data(:,1)) 
hold on 
grid on 
plot(time_, result_x_1 .* -1)
title("Results with F(t) = 100 N")
ylim([-4 4])
legend('x_1 simulation','x_1 lagrange')
subplot(3,1,2)
plot(out.simulation_results.Time(:), out.simulation_results.Data(:,2)) 
hold on 
grid on 
plot(time_, result_x_2)
legend('x_2 simulation','x_2 lagrange')
subplot(3,1,3)
plot(out.simulation_results.Time(:), out.simulation_results.Data(:,3)) 
hold on 
grid on 
plot(time_, result_theta .* -1)
legend('\theta simulation','\theta lagrange')