How to solve Multiple DOF Mass Spring Damper system and find/plot acceleration, velocity & displacement vs time

I know to solve 1 DOF system but i dont know to solve multiple DOF (Matrix will come instead of single value);

Also how to find acceleration?

I am getting absurd/wrong values. Kindly correct the code

Function Code :-

% ODE Function
function dxdt = Ransom(t, x, M, K, C)
m1 = 25.81;m2 = 15.36; m3 = 40.35; m4 = 10.12;   
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620;
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000;
w = 50; F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
dxdt = zeros(8, 1) ;
dxdt(1) = x(5) ;
dxdt(2) = x(6) ;
dxdt(3) = x(7) ;
dxdt(4) = x(8) ;
dxdt(5) = 5*sin(w*t) -K(1)/M(1) * x(1) - K(2)/M(1) * x(2) -C(1)/M(1) * x(5) - C(2)/M(1) * x(6) ;
dxdt(6) = -K(5)/M(6) * x(1) - K(6)/M(6) * x(2) - K(7)/M(6) * x(3) -C(5)/M(6) * x(5) - C(6)/M(6) * x(6) - C(7)/M(6) * x(7) ;
dxdt(7) = -K(10)/M(11) * x(1) - K(11)/M(11) * x(2) - K(12)/M(11) * x(3) -C(10)/M(11) * x(5) - C(11)/M(11) * x(6) - C(12)/M(11) * x(7) ;
dxdt(8) = -K(15)/M(16) * x(3) - K(16)/M(16) * x(4) -C(15)/M(16) * x(7) - C(16)/M(16) * x(8) ;
end

Calling Code: -

clear all;clc;
m1 = 25.81;m2 = 15.36; m3 = 40.35; m4 = 10.12;   
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620;
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000;
w = 50; F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
tspan = [0 1000] ;
y0 = [0 0 0 0 0 0 0 0] ;
[t,x]=ode45('Ransom',tspan,y0);
x_ = x(:, 1:4) ;
xdot_ = x(:, 5:8) ;
subplot(2,1,1)
plot(t, x_);
subplot(2,1,2)
plot(t, xdot_);

Answer 1

[AJMIT KUMAR]

Finally I did it myself

Function Code: -

% ODE Function
function dxdt = Ransom(t,x)
m1 = 25.81; m2 = 15.36; m3 = 40.35; m4 = 10.12;       % Mass(kg)
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620;            % Damping Coefficient 
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000;   % Stiffness
w = 100;                                         % Angular Velocity (rad/s) 
F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
dxdt = zeros(8,1) ;
% x(1),x(2),x(3)and x(4) are displacement of mass m1, m2, m3 & m4 respectively
dxdt(1) = x(5) ;      % velocity of mass m1
dxdt(2) = x(6) ;      % velocity of mass m2
dxdt(3) = x(7) ;      % velocity of mass m3
dxdt(4) = x(8) ;      % velocity of mass m4
dxdt(5) = (F*sin(w*t))/M(1) -K(1)/M(1) * x(1) - K(2)/M(1) * x(2) -C(1)/M(1) * x(5) - C(2)/M(1) * x(6) ;
dxdt(6) = -K(5)/M(6) * x(1) - K(6)/M(6) * x(2) - K(7)/M(6) * x(3) -C(5)/M(6) * x(5) - C(6)/M(6) * x(6) - C(7)/M(6) * x(7) ;
dxdt(7) = -K(10)/M(11) * x(2) - K(11)/M(11) * x(3) - K(12)/M(11) * x(4) -C(10)/M(11) * x(6) - C(11)/M(11) * x(7) - C(12)/M(11) * x(8) ;
dxdt(8) = -K(15)/M(16) * x(3) - K(16)/M(16) * x(4) -C(15)/M(16) * x(7) - C(16)/M(16) * x(8) ;
end

Calling Code: -

clear all;clc;
% Displacement and Velocity Calculation
tspan = [0 1];             % Time (sec)
y0 = [1 1 1 1 0 0 0 0];    % Initial conditions
[t,x]=ode45('Ransom',tspan,y0);
x_ = x(:,1:4);             % All Displacement
xdot_ = x(:,5:8);          % All Velocity

% Plot Displacement wrt Time
subplot(3,1,1)
plot(t, x_ ,'color','b','LineWidth',2);
grid on
xlabel('Time (sec)')
ylabel('Displacement(m)')
title('Displacement Vs Time')

% Plot Velocity wrt Time
subplot(3,1,2)
plot(t, xdot_ ,'color','g','LineWidth',2);
grid on 
xlabel('Time (sec)')
ylabel('Velocity (m/s)')
title('Velocity Vs Time')

% Acceleration Calculation
m1 = 25.81; m2 = 15.36; m3 = 40.35; m4 = 10.12;       % Mass(kg)
c1 = 3564; c2 = 5685; c3 = 8550; c4 = 620;            % Damping Coefficient 
k1 = 950000; k2 = 262800; k3 = 213000; k4 = 450000;   % Stiffness
w = 100;                                         % Angular Velocity (rad/s) 
F = 5;
M = [m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4];
C = [c1+c2 -c2 0 0; -c2 c2+c3 -c3 0; 0 -c3 c3+c4 -c4; 0 0 -c4 c4];
K = [k1+k2 -k2 0 0; -k2 k2+k3 -k3 0; 0 -k3 k3+k4 -k4; 0 0 -k4 k4];
a1 = (F*sin(w*t))/M(1) -K(1)/M(1) * x(:,1) - K(2)/M(1) * x(:,2) -C(1)/M(1) * x(:,5) - C(2)/M(1) * x(:,6) ;
a2 = -K(5)/M(6) * x(:,1) - K(6)/M(6) * x(:,2) - K(7)/M(6) * x(:,3) -C(5)/M(6) * x(:,5) - C(6)/M(6) * x(:,6) - C(7)/M(6) * x(:,7) ;
a3 = -K(10)/M(11) * x(:,2) - K(11)/M(11) * x(:,3) - K(12)/M(11) * x(:,4) -C(10)/M(11) * x(:,6) - C(11)/M(11) * x(:,7) - C(12)/M(11) * x(:,8) ;
a4 = -K(15)/M(16) * x(:,3) - K(16)/M(16) * x(:,4) -C(15)/M(16) * x(:,7) - C(16)/M(16) * x(:,8) ;
a = [a1 a2 a3 a4];

% Plot Acceleration wrt Time
subplot(3,1,3)
plot(t, a ,'color','r','LineWidth',2);
grid on 
xlabel('Time (sec)')
ylabel('Acceleration (m/s^2)')
title('Acceleration Vs Time')

Answer 2

[Stephan]

tspan=[0 10];   % Simulation time
z0=[0;0];       % Initial conditions for forced vibration
[t,z]=ode45(@forced_vibration,tspan,z0);
subplot(2,1,1)
plot(t,z(:,1),'color','b','LineWidth',2);
grid on
xlabel('Time (t) (sec)')
ylabel('Displacement (x) (m)')
title('Displacement Vs Time (Forced Vibration)')
subplot(2,1,2)
plot(t,z(:,2),'color','b','LineWidth',2);
grid on 
xlabel('Time (t) (sec)')
ylabel('Velocity (xdot) (m/s)')
title('Velocity Vs Time (Forced Vibration)')

function f = forced_vibration(t,z)
m = 10;         % Mass(kg)
c = 10;         % Damping coefficient (Ns/m)
k = 1000;       % Stiffness coefficient (N/m)
F = 1;          % Force (N)
w = 50;         % Excitation frequency(rad/s)
f = [z(2);((F/m)*sin(w*t)-(c/m)*z(2)-(k/m)*z(1))];
end

[AJMIT KUMAR]

@Stephan Beg your pardon I didnt provide the value/data for Multiple DOF, now provided, kindly solve for Multiple DOF, the example i have given is for single DOF that I know

[AJMIT KUMAR]

@Stephan Now I have written the code but I am getting wrong values, kindly look on it