sd_id=hdfsd('start','./out/q00000.hdf','read');
sds_index=0;
sds_id=hdfsd('select',sd_id,sds_index);
start=[0 0]; stride=[1 1]; edge=[32 32];
[q0,status]=hdfsd('readdata',sds_id,start,stride,edge);
status = hdfsd('end',sd_id);

sd_id=hdfsd('start','./out/q00001.hdf','read');
sds_index=0;
sds_id=hdfsd('select',sd_id,sds_index);
start=[0 0]; stride=[1 1]; edge=[32 32];
[q1,status]=hdfsd('readdata',sds_id,start,stride,edge);
status = hdfsd('end',sd_id);

mx=32; my=32;
dt=0.5; dx=1; dy=1; dtdx=dt/dx; dtdy=dt/dy;
qLF=zeros([mx my 21]);
qLF(:,:,1)=q0;
u=1; v=0.5; % Matlab transposes x,y

for n=2:21
  contourf(qLF(:,:,n-1))
  axis('square')
  pause(0.5)
  for i=2:mx-1
      for j=2:my-1
          % Lax-Friedrichs
          %f1=0.5*(u*qLF(i,j,n-1)+u*qLF(i+1,j,n-1))-0.25*dx/dt*(qLF(i+1,j,n-1)-qLF(i,j,n-1));
          %f0=0.5*(u*qLF(i-1,j,n-1)+u*qLF(i,j,n-1))-0.25*dx/dt*(qLF(i,j,n-1)-qLF(i-1,j,n-1));
          %g1=0.5*(v*qLF(i,j,n-1)+v*qLF(i,j+1,n-1))-0.25*dy/dt*(qLF(i,j+1,n-1)-qLF(i,j,n-1));
          %g0=0.5*(v*qLF(i,j-1,n-1)+v*qLF(i,j,n-1))-0.25*dy/dt*(qLF(i,j,n-1)-qLF(i,j-1,n-1));          
          % Lax-Wendroff
          f1=0.5*(u*qLF(i,j,n-1)+u*qLF(i+1,j,n-1))-0.5*dtdx*u*(qLF(i+1,j,n-1)-qLF(i,j,n-1)); %-0.125*dtdy*v*(qLF(i,j+1,n-1)-qLF(i,j-1,n-1)+qLF(i+1,j+1,n-1)-qLF(i+1,j-1,n-1));
          f0=0.5*(u*qLF(i-1,j,n-1)+u*qLF(i,j,n-1))-0.5*dtdx*u*(qLF(i,j,n-1)-qLF(i-1,j,n-1)); %-0.125*dtdy*v*(qLF(i-1,j+1,n-1)-qLF(i-1,j-1,n-1)+qLF(i,j+1,n-1)-qLF(i,j-1,n-1));          
          g1=0.5*(v*qLF(i,j,n-1)+v*qLF(i,j+1,n-1))-0.5*dtdy*v*(qLF(i,j+1,n-1)-qLF(i,j,n-1)); %-0.125*dtdx*u*(qLF(i+1,j,n-1)-qLF(i-1,j,n-1)+qLF(i+1,j+1,n-1)-qLF(i-1,j+1,n-1));
          g0=0.5*(v*qLF(i,j-1,n-1)+v*qLF(i,j,n-1))-0.5*dtdy*v*(qLF(i,j,n-1)-qLF(i,j-1,n-1)); %-0.125*dtdx*u*(qLF(i+1,j-1,n-1)-qLF(i-1,j-1,n-1)+qLF(i+1,j,n-1)-qLF(i-1,j,n-1));
          qLF(i,j,n)=qLF(i,j,n-1)-dtdx*(f1-f0)-dtdy*(g1-g0);
      end
  end  
end