If the system were bounded-input-bounded-output (BIBO) stable, then the steady state output in response to input y(t) = A*sin(w*t) would be zss(t) = M*A*sin(wt + phi), where M and phi are determined by the magnitude and phase of the system transfer function evaluated at s = 1j*w.

However, this system has no damping and is not BIBO stable, as indicated by the eigenvalues of the A-matrix all on the imaginary axis

A = [0,1,0,0;-64,0,64,0;0,0,0,1;64,0,-128,0];

e = eig(A)

e =

0.000000000000000e+00 + 1.294427190999917e+01i
0.000000000000000e+00 - 1.294427190999917e+01i
0.000000000000000e+00 + 4.944271909999160e+00i
0.000000000000000e+00 - 4.944271909999160e+00i

So, in this case we have to consider two cases:

- w not equal to 12.94427.... and w not equal to 4.94427.... In this case, the output of the system will be the sum of zss as defined above and and two other sine waves at frequencies 12.944... an 4.94427.... The amplitude and phase of these two other sine waves would have to be determined based on A and w.
- w = 12.9442.... or w = 4.94427.... In this case the output will include a term like C*t*sin(w*t), i.e., the output will grow indefinitely

Finding closed form expressions in either case is feasible, but could be painful.

In either case, you could use lsim to approximate the output via simulation; make sure the the time vector has a time separation small relative to the larger of 12.944 or w, mabye something like dt < 1/10/max(w,12.9444) As expected, the output is the sum of three sinusoids at the frequencies determined by w and the imaginary part of e.

w = (0:(numel(Y)-1))/numel(Y)*2*pi/dt;

w = abs(imag(e(1))); a = 1;

This output is dominated by the growing sine wave at 12.444 rad/sec, but it does have a 4.944 rad/sec component

w = (0:(numel(Y)-1))/numel(Y)*2*pi/dt;

This analysis would be much simpler if you added a dashpot (or damper) in parallel with two springs.