As your ODE is written above we have to ask: Why bother? y is deterministically given provided your input parameters and the time. It should be trivial to calculate it for the times in you tn output. In case you have a slightly more complicated driving function you could perhaps do something like this:
function dz = system(t,z,option,K,M,fe,y_fcn)
dz=zeros(2,1);
y = y_fcn(t); %base excitation
dz(1)=z(2);
dz(2)=(K*(y-z(1))/M;
end
This would let you send in "any" function to your ODE:
% Linear Parameters________________________________________________________
K=5000;% linear stiffness
M=1; % system's mass
%Input Parameters__________________________________________________________
y0 = 1; % amplification of base input
T = 10;
fe = [1 2 3 4];
n_freq = length(fe);
z0 = 0; % you'll have to set this
for kk=1:n_freq
f_current = @(t) y0*sin(2*pi*fe(kk)*t);
[tn,z] = ode15s(@(t,z)system(t,z,K,M,fe(kk),f_current),[t(1) t(end)],z0);
y_exc{kk} = f_current(tn);
end
HTH
