% Wave packet for the Lax-Wendroff eq.
x=[-1:0.05:1]; x1=[-2:0.05:-1]; x2=[1:0.05:9];

y=cos(5*pi*x).*(cos(pi*x/2)).^2;

yy=[x1*0 y x2*0];
xx=[x1 x x2];

clf;
plot(xx,yy);

hold on;
plot(x,(cos(pi*x/2)).^2,'-.');
plot(x,-(cos(pi*x/2)).^2,'-.');

pause



h=0.05;
z1=yy;
z2=z1;
z=z1;

lambda=0.95; k=lambda*h;

% z2=[x1*0, cos(5*pi*(x-k)).*(cos(pi*(x-k)/2)).^2, x2*0];


%for i=2:length(z1)-1
%   z2(i)=(z1(i+1)+z1(i-1))/2-lambda*(z1(i+1)-z1(i-1))/2;
%end

%plot(xx,z,'r');
%pause
   
t=0;
while t*lambda*h<8
   t=t+1;
   for i=2:length(z1)-1
      z(i)=z1(i)-lambda/2*(z1(i+1)-z1(i-1))+lambda^2/2*(z1(i+1)-2*z1(i)+z1(i-1));
   end
   z1=z;
   z2=z;
   
   clf
   plot(xx,yy);

   hold on;
   plot(x,(cos(pi*x/2)).^2,'-.');
   plot(x,-(cos(pi*x/2)).^2,'-.');
   
%   axis([-2 10 -3 3]);
   
   plot(xx,z,'r');
   
   
   % true packet moving
   tt=t*k; r=-1+tt:0.05:1+tt; plot(r,(cos(pi*(r-tt)/2)).^2,'-.');
   
   
   % phase speed
   alpha=1-(h*5*pi)^2/6*(1-lambda^2);
   tt=alpha*t*k; r=-1+tt:0.05:1+tt; 
%   plot(r,cos(5*pi*(r-tt)).*(cos(pi*(r-tt)/2)).^2,'g');
   
   % computed packet moving
   % group speed 
   gamma=alpha-1/3*(h*5*pi)^2*(1-lambda^2);
   tt=gamma*t*k; r=-1+tt:0.05:1+tt; plot(r,(cos(pi*(r-tt)/2)).^2,'r-.');
   
   tt2=alpha*t*k; r=-1+tt:0.05:1+tt; 

   plot(r,cos(5*pi*(r-tt2)).*(cos(pi*(r-tt)/2)).^2,'g');

   
   pause

end