N=2;

w0=.3;

w0sq=w0^2;
wsq=w0sq*repmat(1,[N+1 1]);

A=coup_osc_oper(wsq,0);

[V,D]=eigs(A,N); % normal modes

% eigen frequencies
sig=sqrt(-diag(D));

