NDOF Linear Shear Building - Forcing Function problem set up

I have developed this example where a simple NDOF shear building is analyzed under Step Load [P], for damped/undamped condition using newmark's method. However, I am having difficulties setting up the Forcing Matrix as F(x)=F1+F2*Sin(Omega*t); Right now it is just a step load matrix, with pre-defined forces in each floor.

        if true
          % code 
        end  

% Code by :        Azizi                                                  |
%                  CEE 536 - Structural Dynamics                          |
%                                                                         |
%--------------------------------------------------------------------------
%
% This MATLAB program solves the NDOF linear shear building model using 
% newmark's numerical integration method. In order to use the program;
% first describe the initial conditions and bounder's as matrixes. Values
% to input are ndof,M,K,C,P,and type of newmark acceleration. 
% 
%--------------------------------------------------------------------------
clear ;clc ;

EI = 5*10^9 ;
h = 100 ;
K = (EI/h^3)*[6     -4    0    0    0 ;       % K - Stiffness Matrix (INPUT)
             -4     10   -6    0   0 ;       
              0     -6   11   -4   0 ;
              0      0   -4   6  -2 ; 
              0      0    0   -2   2];

   ndof = length(K) ;         
   m = 200 ;
   M = m*diag([2 1 1 1 1]) ;                 % M - Mass Matrix  (INPUT)

   %%%%%%%Modal Analysis                     
  [V,D]=eig(K,M);
  [W,k]=sort(diag(D));
  V=V(:,k); 
  Factor=diag(V'*M*V);
  Phi=V*inv(sqrt(diag(Factor)));
  Omega=diag(sqrt(Phi'*K*Phi)); 
  w = diag(Omega(1:ndof)) ;

%-------------------------------------------------------%
%%%%%%%%%%%%%Provide damping option &&&&&&&

sys = 'damped' ;       %Select if system Damped
switch sys
    case 'undamped'     %  C - If System undamped
        C = zeros(size(K)) ;
    case 'damped'
        C= 2*w*0.2;     % C - Damping Matrix (Damped)
end

%-------------------------------------------------------
%%%%%%%%%%%%%NEWMARK'S METHOD : LINEAR SYSTEM &&&&&&&
acc = 'Linear' ;   % acc - Type of Newmark's Method to be used 
switch acc
    case 'Average' 
        gama = 1/2 ;beta = 1/4 ;
    case 'Linear'
        gama = 1/2 ;beta = 1/6 ;
end
% Time step
ti = 0. ;
tf = 20. ;
dt = 0.1 ;
t = ti:dt:tf ;
nt = fix((tf-ti)/dt) ;
n = length(M) ;

% Constants used in Newmark's integration
a1 = gama/(beta*dt) ; a2 = 1/(beta*dt^2) ;
a3 = 1/(beta*dt) ;       a4 = gama/beta ;
a5 = 1/(2*beta) ;        a6 = (gama/(2*beta)-1)*dt ;
P = 100*[10 20 30 40 50]' ;   % P - Force MAtrix

disp = zeros(n,nt) ;
vel = zeros(n,nt) ;
accl = zeros(n,nt) ;
U = zeros(ndof,nt) ;

% Initial Conditions for displacement and acceleration
% Individual input for initial position of each building can be made:
disp(:,1) = zeros ; %Building is assumed at rest
vel(:,1) = zeros ;
U(:,1) = Phi*disp(:,1) ;
P0 = zeros(n,1) ;
accl(:,1) = M\(P-C*vel(:,1)-K*disp(:,1)) ;

Kcap = K+a1*C+a2*M ;

a = a3*M+a4*C ;
b = a5*M+a6*C ;

% Tme step starts
for i = 1:nt
    delP = P0+a*vel(:,i)+b*accl(:,i) ;
    delu = Kcap\delP ;
    delv = a1*delu-a4*vel(:,i)-a6*accl(:,i) ;
    dela = a2*delu-a3*vel(:,i)-a5*accl(:,i);
    disp(:,i+1) = disp(:,i)+delu ;      % disp - modal displacement's 
    vel(:,i+1) = vel(:,i)+delv ;        % vel - modal velocities
    accl(:,i+1) = accl(:,i)+dela ;      % accl - modal accelerations
    U(:,i+1) = Phi*disp(:,i+1) ;        % U - system's displacement

end

%%%%%%%%GRAPHICAL PLOT%%%%%%
figure(1);
subplot(2,3,1);
plot(t,disp(1,:),'-','Color','blue') ;
xlabel('Time in sec') ; ylabel('displacement');
title('First Floor Displacement');

subplot(2,3,2);
plot(t,disp(2,:),'-','Color','red');
xlabel('Time in sec') ; ylabel('displacement');
title('Second Floor Displacement');

subplot(2,3,3)
plot(t,disp(3,:),'-','Color','yellow');
xlabel('Time in sec') ; ylabel('displacement');
title('Third Floor Displacement');

subplot(2,3,4)
plot(t,disp(4,:),'-','Color','green');
xlabel('Time in sec') ; ylabel('displacement');
title('Fourth Floor Displacement');

subplot(2,3,5)
plot(t,U(5,:),'-','Color','black');
xlabel('Time in sec') ; ylabel('displacement');
title('Building Displacement');