% Wave packet for Leap-frog
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,'-.');





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

% Initialize Leap frog using Lax-Friedrichs
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*(z2(i+1)-z2(i-1));
   end
   z1=z2;
   z2=z;
   
   clf
   plot(xx,yy);

   hold on;
   plot(x,(cos(pi*x/2)).^2,'-.');
   plot(x,-(cos(pi*x/2)).^2,'-.');
   
   
   % numerical solution
   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,'-.');
   
   % computed packet moving
   % group speed = 0.9545
   tt=0.9545*t*k; r=-1+tt:0.05:1+tt; plot(r,(cos(pi*(r-tt)/2)).^2,'r-.');
   
   % phase speed = 0.9872
   tt2=0.9872*t*k; r=-1+tt:0.05:1+tt; 
%   tt=0.9545*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



plot(xx,z,'r');