%% Example MDOF calls
% Based on example in Chopra, 3rd Edn., p 458
clear all;

% Define system matrices
M = 1/386.*diag( [400, 400, 200] );
C = [0.98304, 0.32069, 0.27191;
     0.32069, 1.52690, 0.32069;
     0.27191, 0.32069, 0.62748];
K = 610 .* [2,  -1, 0;
            -1, 2,  -1;
            0,  -1, 1;];
nDOF = 3;

% Create time vector
t0 = 0;
dt = 0.01;
tf = 5;
t = t0:dt:tf;

% Initial conditions
u0 = zeros(nDOF,1);
du0 = zeros(nDOF,1);

% Define forcing function
amp = 5000;
omega = 12;
n = length(t);
F = zeros(n,nDOF);
F(:,1) = amp*cos(omega*t);  % load at 1st DOF

[u du ddu] = Newmark_Beta(t, M, C, K, F, u0, du0);

figure;
subplot(2,2,1:2);   
    plot(t,u);    
    xlabel('Time (s)'); 
    ylabel('Displacement');
subplot(2,2,3);     
    plot(t,du); 
    xlabel('Time (s)'); 
    ylabel('Velocity');
subplot(2,2,4);     
    plot(t,ddu);            
    xlabel('Time (s)'); 
    ylabel('Acceleration');