function wavepacket
clear; close all;
format long;
%Parameters
Ao=1;
a=2.1e6;
pw=.005;
%Time window in units of laser periods.
tmin=0;
tmax=50;
%number of frames
nframe=400;
dt=(tmax-tmin)/nframe;
%number of pixels in a wavelength
nmax=100;
D=1/nmax;
%Define frame
x=-1:D:1;
z=0:D:10;
[X,Z]=ndgrid(x,z);
for n=1:nframe+1;
  t=tmin+(n-1)*dt;
  rho=Evaluate(X,Z,t,pw,Ao,a);
  peak=max(max(rho));
  rho=rho/peak;
  colormap([.25 0 .5;.25 0 .75;.25 0 1;0 0 1;0 .25 1;0 .5 1;0 .75 1;0 .9 .9;
      0 .9 .7;0 .95 .3;0 .97 0;.3 1 0;.7 1 0;.87 1 0; 1 1 0;1 .87 0;1 .75 0;
      1 .5 0;1 .25 0;1 0 0;1 .25 .25;1 .5 .5;1 .75 .75;1 1 1]);
  shading interp
  imagesc(rho)
  set(0,'DefaultImageCreateFcn','axis image')
  drawnow
  M(n)=getframe;
end;
movie2avi(M,'volkov')

function rho=Evaluate(X,Z,t,pw,Ao,a);
pw2=pw^2;
pw3=pw^3;
pw4=pw^4;
Ao2=Ao^2;
pw4at=pw4*a*t;
fxpw3=pw3*a*Ao/(2*pi);
fpw3=pw^3*a*Ao2*t/4;
fprintf('t = %g    pw^4*a*t = %g     pw^3*f = %g     pw^3*fx = %g\n',t,pw4at,fpw3,fxpw3)
fx=a*Ao/(2*pi)*sin(2*pi*(Z-t));
fxd=-a*Ao*cos(2*pi*(Z-t));
f=a*Ao2*((Z-t)/4+sin(4*pi*(Z-t))/(16*pi));
fd=-a*Ao2*(cos(2*pi*(Z-t))).^2/2;
ntp=a*t+f;
ntpd=a+fd;
ntm=a*t-f;
ntmd=a-fd;
nx=a*X+fx;
nxd=fxd;
nz=a*Z+f;
nzd=fd;
al=1+i*pw2*ntp;
ald=i*pw2*ntpd;
b=1+i*pw2*ntm+pw4*fx.^2./al;
bd=i*pw2*ntmd+2*pw4*fx.*fxd./al-pw4*ald.*(fx./al).^2;
n=nz+i*pw2*fx.*nx./al;
nd=nzd-i*pw2*ald.*fx.*nx./(al.^2)+i*pw2*(fxd.*nx+fx.*nxd)./al;
q=real(nx.^2./al+n.^2./b);
w=imag((1-pw2*nx.^2./(2*al)).*ald./al+(1-pw2*n.^2./b).*bd./(2*b)+pw2*(nx.*nxd./al+n.*nd./b));
rho=(a*pw/sqrt(pi))^3/a*exp(-pw2*q).*(ntmd+w)./((abs(al)).^2.*abs(b));