function [V] = vfunc(Y,x,xdot,omegaN,alpha)

V=[0,1;-omegaN^2,-2*alpha]*Y+[0;1]*(omegaN^2*x+2*alpha*xdot);