%Matlab Code for Two Degree of freedom vibratory System by Using 4th order Runge Kutta

function RK_Two_masses
global t dt
tf = 20;  %(sec)
dt =0.001; %(sec) 

%initial conditions
x_0_1=0; % Initial displacement of the first mass
x_0_1_dot=0; % Initial velocity of the first mass
x_0_2=0; % Initial displacement of the second mass
x_0_1_dot=0; % Initial velocity of the second mass

xstate= [x_0_1; x_0_1_dot; x_0_2; x_0_1_dot];

%calculations

for i=1:tf/dt
t = dt*i;    
xstate = RK(xstate);

%results
z1(i) = xstate(1); % (m)
z2(i) = xstate(2); % (m)
end
figure(1);plot(0:dt:tf,[0,z1*1000],'k-',0:dt:tf,[0,z2*1000],'r:');
xlabel('time (sec)'); ylabel('response (mm)');
legend('X_1 (t)','X_2 (t)');


function xstate= RK(xstate)
global dt
f1 = dt * x_dot( xstate       );
f2 = dt * x_dot( xstate + f1/2);
f3 = dt * x_dot( xstate + f2/2);
f4 = dt * x_dot( xstate + f3  );
xstate = xstate + 1/6 *(f1 + 2*f2 + 2*f3 + f4);


function xdot=x_dot(X)

% SOME ARBITRARY coefficients 
 % (no damping)


M = [9, 0; 0, 1];
K = [27,-3; -3, 5];

% define F1 and zv
F = [F1;0];

A =[zeros(2,2), eye(2); -inv(M)*K, zeros(2,2)];
b =[0;0;inv(M)*F];

xdot = A*X+b;

function F1_t = F1
global t
F1_t = 10 * sin(10*t);
