n=256;
r=8;
D = .5;
tmax = 2*n;
u0=u=v=zeros(n,1);
for i=1:n
	u0(i) = exp(-(i-n/2)^2 * .001)*sin(i*.8);
end
u = v = u0;

%Laplace
L = toeplitz([-2.0 1 zeros(1,n-2)]);
L(1,n)=1;
L(n,1)=1;

S = [];
for t=0:tmax
	w = -v + 2*u + D*L*u;
	v = u;
	u = w;
	
	% draw the animation
	clf;
	plot(u);
	axis([0 n -.5 1.5]);
	drawnow;
	if mod(t,5) == 0
		S = [S u];
	end
end

[U,Sigma,VT] = svd(S);
%plot(U(:,1:3)) % Some of the dominant modes

% now the ROM waves
Phi = U(:,1:r);
r0 = u0;
ru0 = Phi' * r0;
ru = rv = ru0;
rL = Phi'*L*Phi;

u=v=u0;

for t=1:tmax
	rw = -rv + 2*ru + D*rL*ru;
	rv = ru;
	ru = rw;
	
	w = -v + 2*u + D*L*u;
	v = u;
	u = w;

	clf;
	p1 = plot(Phi*ru); 
	set(p1,'color','red');

	hold on;

	p2 = plot(u);
	axis([0 n -.2 1.1]);
	set(p1,'linewidth',2);
	set(p2,'linewidth',2);
	drawnow;

	outpath = sprintf('output/out.%03d.png',t);
	print(outpath);
end