Problem with Cavity Driven Flow using Lattice Boltzmann Method

I have written a code for LBM D2Q9 lattice structure, but I have some problems with convergence and pressure contours

if true
  %%2D Lid Driven Cavity Flow %%
nx = 100; ny = 100;                                  % Mesh Size %
maxT = 5e3;                                          % Time Steps %
u_ini = 0.1; v_ini = 0;                             % Initial Velocities %
Re = 0100;                                           % Reynolds Number %
% tou = (1/2)*(((6*u_ini) / (Re)) + 1);                % Parameters % 
tou=1;
%%D2Q9 Lattice Constants %%
n = 9;                                               % Lattice Points %
w1 = 4/9;                                            % Direction - Origin %
w2 = 1/9;                                            % Direction - CH, CV %
w3 = 1/36;                                           % Direction - Diagonal %
rho = 2.7;                                             % Initial Density %
cs2 = 1/3;                                           % Constant %
X = [w2 w3 w2 w3 w2 w3 w2 w3 w1];
f = reshape( X'*ones(1,nx*ny),nx,ny,n);              % Initial Equilibrium %
f_eq = f;
for t = 1:maxT
%%Streaming Process %%
  f(:,:,1) = f([nx 1:nx-1],:,1); 
  f(:,:,2) = f([nx 1:nx-1],[ny 1:ny-1],2); 
  f(:,:,3) = f(:,[ny 1:ny-1],3);
  f(:,:,4) = f([2:nx 1],[ny 1:ny-1],4);
  f(:,:,5) = f([2:nx 1],:,5);
  f(:,:,6) = f([2:nx 1],[2:ny 1],6);
  f(:,:,7) = f(:,[2:ny 1],7); 
  f(:,:,8) = f([nx 1:nx-1],[2:ny 1],8);
%%Boundry Conditions %%
  f(1,:,1) = f(1,:,5); 
  f(1,:,2) = f(1,:,6); 
  f(1,:,8) = f(1,:,4);
  f(nx,:,4) = f(nx,:,8); 
  f(nx,:,5) = f(nx,:,1); 
  f(nx,:,6) = f(nx,:,2);
  f(:,1,2) = f(:,1,6); 
  f(:,1,3) = f(:,1,7); 
  f(:,1,4) = f(:,1,8);
  rhon = f(:,ny,9)+f(:,ny,1)+f(:,ny,5)+2*(f(:,ny,3)+f(:,ny,4)+f(:,ny,2));
  f(:,ny,7) = f(:,ny,3);
  f(:,ny,6) = f(:,ny,2)+0.5*(f(:,ny,1)-f(:,ny,5))-u_ini.*rhon/2;
  f(:,ny,8) = f(:,ny,4)+0.5*(f(:,ny,5)-f(:,ny,1))+u_ini.*rhon/2;
%%Macroscopic Variables Computation %%
  rho = sum(f,3);
  u = (sum(f(:,:,[1 2 8]),3) - sum(f(:,:,[4 5 6]),3))./rho;
  v = (sum(f(:,:,[2 3 4]),3) - sum(f(:,:,[6 7 8]),3))./rho;
  P = rho*cs2;
%%New Equilibrium Computation %%
  f_eq(:,:,1) = w2*rho.*(1+3*u+9/2*u.^2-3/2*(u.^2+v.^2));
  f_eq(:,:,3) = w2*rho.*(1+3*v+9/2*v.^2-3/2*(u.^2+v.^2));
  f_eq(:,:,5) = w2*rho.*(1-3*u+9/2*u.^2-3/2*(u.^2+v.^2));
  f_eq(:,:,7) = w2*rho.*(1-3*v+9/2*v.^2-3/2*(u.^2+v.^2));
  f_eq(:,:,2) = w3*rho.*(1+3*(u+v)+9/2*(u+v).^2-3/2*(u.^2+v.^2));
  f_eq(:,:,4) = w3*rho.*(1+3*(-u+v)+9/2*(-u+v).^2-3/2*(u.^2+v.^2));
  f_eq(:,:,6) = w3*rho.*(1-3*(u+v)+9/2*(u+v).^2-3/2*(u.^2+v.^2));
  f_eq(:,:,8) = w3*rho.*(1+3*(u-v)+9/2*(u-v).^2-3/2*(u.^2+v.^2));
  f_eq(:,:,9) = w1*rho.*(1-3/2*(u.^2+v.^2));
%%Collision%%
  f = ((1-tou)*f) + (tou*f_eq);   
%%Normalizing Velocity Scale %%
  u1=u/u_ini;
  v1=v/u_ini;
  if t > 1
      bb = (u1-a1).*(u1-a1);
      cc = sum(sum(bb));
      dd = (v1-a2).*(v1-a2);
      ee = sum(sum(dd));
      Erru = abs(sqrt(cc + ee));
      u12 = u1.*u1;
      v12 = v1.*v1;
      ff = sum(sum(u12));
      gg = sum(sum(v12));       
      Erruv = abs(sqrt(ff + gg));
      Err(t,1) = Erru/Erruv;
      a1 = u1;
      a2 = v1;
  else
      a1 = u1;
      a2 = v1;
      Err(t,1) = 0;                                                            %#ok<*SAGROW>
  end
end
%%Plotting Figures %%
figure (1)
contour(u1',nx);
title('Streamlines')
xlabel('nx');
ylabel('ny');
figure(2)
plot(1:maxT,Err);
title('Error Curve')
xlabel('Iterations');
ylabel('Error');
figure(3)
plot(u1(nx/2,:),0:1/nx:(1-(1/nx)));
title('Horizontal Velocity Profile')
xlabel('X');
ylabel('Y');
figure(4)
plot(0:1/ny:(1-(1/ny)),v1(:,ny/2));
title('Vertical Velocity Profile')
xlabel('X');
ylabel('Y');
figure(5)
contour(P',nx);
title('Pressure Contour')
end