Multi-degree Forced spring-mass system with damper energy conservation

Question

Edit: It was kindly pointed out to me in the comments that because of the MATLAB code I should post in superuser.com as this is a mathematics-oriented forum. However I want to clarify that there is no problem with the code, its just the theory behind it that I am not comprehending. The code is there to help with visualisng the problem and getting familiar with the details.

The Problem:

Having writen a MATLAB code to solve a 3 Degree of freedom, forced, spring-mass system with a damper, I try to evaluate the energy balance for every time $t$ of the solution. I am lead to believe from my results that I have done this correctly as the energy balance is virtually $0$ with an input force phase angle of $\phi=[120^{\circ}, 240^{\circ}, 360^{\circ}]^T$. However, by changing the input force phase angles to something arbitraty such as $\phi=[52^{\circ}, 114^{\circ}, 171^{\circ}]^T$, the energy balance does not hold anymore which shouldn't be the case obviously.

Some Context:

The system is modelled as shown in this image, with the following equation:

$\mathbf{M\ddot{x} + C\dot{x} + Kx = F(t)}$

Wehere $\mathbf{F(t)}$ is a simple sinusoidal response. Naturally, the energy in the system is as follows:

$\begin{align} PE & = \frac{1}{2} \mathbf{ x(t)^T K x(t)} && \,\text{i.e. the Potential Energy}\\ KE & = \frac{1}{2} \mathbf{ \dot{x}(t)^T M \dot{x}(t)} && \,\text{i.e. the Kinetic Energy} \\ DE & = \int^t_0 \mathbf{\dot{x}(t)^T C \dot{x}(t)}dt && \,\text{i.e. the Dissipated Energy} \\ EE & = \int^t_0 \mathbf{\dot{x}(t)F(t)}dt && \,\text{i.e. the Excitation Force Energy} \\ \end{align}$

Hence, based on the above, I would assume that the energy balance at any point $t$ would be:

$KE+PE+DE = EE + c$

where $c$ is a constant depending on the energy present from the intial conditions of the system.

Apparently, there seems that there must be a flaw in my logic above, as by running the code attached, by changing the phase angle of the excitation force $\phi$ (line 24) the above does not hold true anymore.

Thanks for your help in advance,

KMT.

P.S.

The Code

clear ; clc

% User Inputs

% Coefficients
m = 3 ;         % Mass, kg
k = 20 ;        % Spring constant, N/m
c = 2 ;         % Damper constant, Ns/m
f0= 10 ;        % Excitation force amplitude, N
f = 0.5;        % Excitation force frequency, Hz

% Degrees of freedom
N = 3 ;

% Solver inputs
tmax = 30 ;     % Maximum solution time
vi = 0 ;        % Initial velocity of 1st DoF

% Matrices
M = diag(repmat(m, 1, N)) ;
K = tridiag(repmat(k, 1, N+1)) ;
C = tridiag(repmat(c, 1, N+1)) ;
F = diag(repmat(f0, 1, N)) ; 
Phi = linspace(2*pi/N,2*pi,N)' ; % <-- Change to [1 2 3]'

% ODE Solution
tspan = [0 tmax] ;
y0 = [zeros(1, N) vi zeros(1, N-1)] ;

[t, y] = ode45(@(t, y) odefcn(t, y, M, K, C, F, Phi, f, N), tspan, y0) ;

[Ec, En, Ee] = energy(t, y, M, K, C, F, Phi, f, N) ;

% Plotting
fig = figure;
gcf = fig ;
fig.Color = 'w' ;
fig.Position = [400 200 900 500] ;

subplot(2, 1, 1)
    yyaxis right
    plot(t, y(:,1:N))
    ylabel('Displacement (m)') ;

yyaxis left
    plot(t, y(:,N+1:2*N))
    ylabel('Velocity (m/s)') ; xlabel('Time (s)') 
    legend(sprintfc('Mass %d',1:N),'Location','Northoutside','Orientation','Horizontal')
    grid on

subplot(2, 1, 2)
    yyaxis right
    plot(t, [Ec En Ee])
    ylabel('Energy (J)') ;

    yyaxis left
    plot(t, sum([Ec En -Ee], 2))
    legend('Total Energy', 'Kinetic + Potential Energy', 'Dissipated Energy','Input Energy','Location','Northoutside','Orientation','Horizontal') 
    ylabel('Energy (J)') ;  xlabel('Time (s)')
    grid on

% Energy function
function [Ec, En, Ee] = energy(t, ymat, M, K, C, F, Phi, f, N) 

% Indexes
i = N ; 
j = N+1 ; 
k = 2*N ; 

% Extract Variables
y  = ymat(:, 1:i)' ;
dy = ymat(:, j:k)' ;

% Energy/Power Equations
KE = 0.5 * dy' * M * dy ;                       % Kinetic energy
PE = 0.5 * y' * K * y ;                         % Potential energy
DP = dy' * C * dy ;                             % Dissipated power
EP = abs(dy' * F*sin(2*pi*f*t' - Phi)) ;        % Excitation power

Ec = diag(KE + PE) ;                            % Sum of potential & kinetic energy
En = cumtrapz(t,diag(DP)) ;                     % Integrated dissipated power over time
Ee = cumtrapz(t,diag(EP)) ;                     % Integrated excitation power

end

% ODE function
function dydt = odefcn(t, y, M, K, C, F, Phi, f, N)

% Indexes
i = N ;
j = N+1 ;
k = 2*N ;

% Preallocate
dydt = zeros(k, 1) ;

% Equation
dydt(1:i) = y(j:k) ;
dydt(j:k) = M \ (-C*y(j:k) -  K*y(1:i) + F*sin(2*pi*f*t - Phi)) ;

end

% Tri-diagonal matrix function
function X = tridiag(x_c)

 x(:,:,1) = diag(x_c(1:end-1)) ;
 x(:,:,2) = diag(x_c(2:end)) ;
 x(:,:,3) = diag(-x_c(2:end-1), -1);
 x(:,:,4) = diag(-x_c(2:end-1), 1) ;

 X = sum(x,3) ;

end