function [u du ddu] = Newmark_Beta(t, M, C, K, F, u0, du0)
% Newmark-Beta Integration for MDOF systems
%
% INPUT
%       [t]:        Time Vector             [n,1]
%       [M]:        System Mass             [nDOF,nDOF]
%       [C]:        System Damping          [nDOF,nDOF]
%       [K]:        System Stiffness        [nDOF,nDOF]
%       [F]:        Applied Load(s)         [n,nDOF]
%       [u0]:       Initial Position        [nDOF]
%       [udot0]:    Initial Velocity        [nDOF]
%
% OUTPUT
%       [u]:        Displacemente Response   [n,nDOF]
%       [du]:       Velocity                 [n,nDOF]
%       [ddu]:      Acceleration             [n,nDOF]
%
% Based on Cheng - Matrix Analysis of Structural Dynamics (p.329), with an 
% allowance is made for the initial acceleration of the system as per
% Chopra, 3rd Edn. p.614) 

gamma   = 1/2;
beta    = 1/4;

nDOF = length(M);
n = length(t);
h = t(2) - t(1);

% Constant terms and effective stiffness
a0 = 1/(beta*h^2);      a1 = gamma/(beta*h);
a2 = 1/(beta*h);        a3 = 1/(2*beta)-1;
a4 = gamma/beta-1;      a5 = h*(gamma/(2*beta)-1);
a6 = h*(1-gamma);       a7 = gamma*h;

Keff = K + a0*M + a1*C;

% Initial acceleration
ddu0 = M\(F(1,:)' - C*du0 - K*u0); % same as inv(M)*(.)

% Initial conditions and output matrices
u    = zeros(n,nDOF);
du   = zeros(n,nDOF);
ddu  = zeros(n,nDOF);
u(1,:) = u0;
du(1,:) = du0;
ddu(1,:) = ddu0;

%   loop for all time steps
for i = 2:n  
    Feff = F(i,:)' + M*(a0*u(i-1,:)' + a2*du(i-1,:)' + a3*ddu(i-1,:)') ...
                     + C*(a1*u(i-1,:)' + a4*du(i-1,:)' + a5*ddu(i-1,:)');
  
    u(i,:)   = Keff\Feff; % Same as inv(Keff)*Feff
    ddu(i,:) = a0*( u(i,:)-u(i-1,:) ) - a2*du(i-1,:) - a3*ddu(i-1,:);
    du(i,:)  = du(i-1,:)  + a6*ddu(i-1,:) + a7*ddu(i,:);
end