function dX = systemsimple(t, X)
% dX = [dQ dP']

M = size(X,1); P = (M/2);
omega = 1; gamma = 1; m = 1;

dX = zeros(M,1); 
% dQ = (1/m)P
dX(1:P)= (1/m)*X(P+1:end);     

% dP' = -kQ-wP
dX(P+1:end) = -m*(omega^2+gamma^2)*X(1:P) - 2* gamma*X(P+1:end);