N=1000;

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));

% Initial conditions
%psi0=[-.3 .3]'; % position
%psidot0=[0 0]'; % velocity
psi0=V(:,900);
psidot0=repmat(0,[N 1]);
% project ICs onto normal modes
a=V'*psi0;
b=V'*psidot0./sig;

% psi = sum_{j=1:N}[phi_j*(a_j*cos(sig_j*t) + b_j*sin(sig_j*t))]
% phi_j = V(:,j) j-th normal mode
% psi=V*(a.*cos(sig*t)+b.*sin(sig*t));

