Do you mean something along these lines:
tspan = [0 2]; % Or whatever you need
IC = [1 0]; % Set your own initial conditions IC = [u(0) dudt(0)]
[t, U] = ode45(@fn, tspan, IC);
u = U(:,1); dudt = U(:,2);
plot(t,u,'r',t,dudt,'b'),grid
function dUdt = fn(t,U)
m=1000;
k=5000;
Ke = 1000; % Set to your own value
Me = 1000; % Ditto
Z = 1000; % Ditto
omega = 10; % Ditto
c=0.02*2*sqrt(Ke*Me);
% I suspect you want z = Z*sin(omega*t), so that
% dz/dt = Z*omega*cos(omega*t)
dzdt = Z*omega*cos(omega*t);
M = [1 0; 0 4*m];
K = [0 -1; 3*k 30*c];
V = [0; 12*c*dzdt];
dUdt = M^-1*(V - K*U);
end
% Let v = u', then your 2nd order equation becomes two ist order equations
% u' - v = 0
% 4*m*v' + 30*c*v + 3*k*u = 12*c*Z*omega*cos(omega*t)
% [1 0][u'] + [0 -1][u] = [0 ]
% [0 4*m][v'] [3*k 30*c][v] [12*c*Z*omega*cos(omega*t)]
