N=1000;

w0=.3;

w0sq=w0^2;
% add random perturbations
rmax=1.05;
rmin=0.93;
r=rand(N+1,1);
r=(rmax-rmin)*r+rmin; % scale
%r=repmat(1,[N+1 1]);
wsq=w0sq*r;

A=coup_osc_oper(wsq,0);

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

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

% Initial conditions
d=1;
n=[0:N+1]';
z=d*n;
n0=round(N/2);
z0=d*n0;
psi0=0.3*exp(-((z-z0).^2)/50).*sin(n*200*pi/(N+1));
psi0=psi0(2:end-1);
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));

