function todam = todam(q,p,t)

% Scattering for Toda lattice with 5 particles on the line.
% The data are a 5-vector q of five positions at time 0, 
% a 5-vector p of momenta at time 0, and time length.
% On the left are the velocities, on the right momenta
% The movie starts at time - t

     [T,Y] = toda(q,p,-t);
  
     N = length(T);

     [T,Y] = toda(Y(N,[1:5]),Y(N,[6:10]),2*t);

     m1 = 1.2*min(min(Y(:,1)),min(Y(:,5))); 

     M1 = 1.2*max(max(Y(:,1)),max(Y(:,5)));

     m2 = 1.2*min(min(Y(:,10)),min(Y(:,6))); 

     M2 = 1.2*max(max(Y(:,6)),max(Y(:,10)));

		  N = length(T);
figure

		  T = T - t;

     for j=1:N
     subplot(1,2,1)
plot(Y([1:j],[1:5]),T([1:j]),'LineWidth',2)
     xlim([m1,M1])
     ylim([-t,t])
     pause(0.1)
     subplot(1,2,2)
plot(Y([1:j],[6:10]),T([1:j]),'LineWidth',2)
     ylim([-t,t])
     xlim([m2,M2])
     pause(0.1)
end

function [T,Y] = toda(q,p,t);


[T,Y] = ode45(@rigid,[0 t], [q p]);



function dy = rigid(t,y)
     N=5;
dy = zeros(2*N,1);  
dy(1)=y(N+1);
dy(N)=y(2*N);
dy(N+1)=-exp(y(1)-y(2));
dy(2*N)=exp(y(N-1)-y(N));
for j=2:N-1
dy(j) = y(j+N);
dy(j+N)=-exp(y(j)-y(j+1))+exp(y(j-1)-y(j));
end


