Lid Driven Cavity Problem

I have simulated Lid Driven Cavity Flow using Lattice Boltzmann Method. But I'm getting "NaN" values for density as well as for velocity also.

Please help me out to fix this.

My code is as follows.

% D2Q9 Model 

clc
clear variables
close all

M = 100;
N = 100;
xE = 1.0; yE = 1.0;
dx = xE/M; dy = yE/N;
x = 0:dx:xE; y = 0:dy:yE;

% Set Reynolds Number = 1000;
nu = 0.01;                                          % Kinematic Viscosity
u0 = 0.1;                                            % Lattice Velocity

Cs = 1/sqrt(3);
w = [1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36,4/9];      % w(9) represents w0 = w(0);
cx = [1,0,-1,0,1,-1,-1,1,0];                        % X-components of Discrete Particle Velocity
cy = [0,1,0,-1,1,1,-1,-1,0];                        % Y-components of Discrete Particle Velocity

omega = 1/(3*nu + 0.5);                             % Equation (8.12)

% Initial Variables
feq = zeros(M+1,N+1,9);                             % f(9) represents rest particle distribution
f = zeros(M+1,N+1,9);   
rho = ones(M+1,N+1);                               % Solution
u = zeros(M+1,N+1);
v = zeros(M+1,N+1);

tol = 1e-1; err = 10; err1 = 0; count = 0;

% Lid Velocity or Velocity of Top Wall
u(:,N+1) = u0;
tic
 for t = 1:1500
% while err>tol
    % Collision
    [f] = collision(M,N,cx,cy,u,v,omega,f,feq,rho,w);
    
    % Streaming
    [f] = stream(f);
    
    [f] = bcs(M,N,u0,f);
%     
    [rho,u,v]=rhoU(M,N,f);


 end
% Plotting

for j = 1:N+1
    u1(j) = u(M/2,j)/u0;
end
for i = 1:M+1
    v1(i) = v(i,N/2)/u0;
end
for i = 1:M+1
    for j = 1:N+1
        U(i,j) = sqrt(u(i,j)^2 + v(i,j)^2);
    end
end
toc
[a,b] = contourf(u')
clabel(a,b)
colorbar;
xlabel('X')
ylabel('Y')
title('Velocity Distribution for Re = 1000')

xx = 0:0.01:1;
yy = 0:0.01:1;

figure(2)
streamline(x,y,v,u,xx,yy)

% -------------------------------------------------------------------------------------
function [f] = collision(M,N,cx,cy,u,v,omega,f,feq,rho,w)

for i = 1:M+1
    for j = 1:N+1
        t1 = u(i,j)*u(i,j) + v(i,j)*v(i,j);
        
        for k = 1:9
            t2 = cx(k)*u(i,j) + cy(k)*v(i,j);
            feq(i,j,k) = w(k)*rho(i,j)*(1 + 3*t2 + (4.5)*t2*t2 - (1.5)*t1);
            f(i,j,k) = (1-omega)*f(i,j,k) + omega*feq(i,j,k);
        end
    end
end
end
% -------------------------------------------------------------------------------------
function [f]=stream(f)
f(:,:,1)=circshift( squeeze(f(:,:,1)), [+1,+0] );
f(:,:,2)=circshift( squeeze(f(:,:,2)), [+0,+1] );
f(:,:,3)=circshift( squeeze(f(:,:,3)), [-1,+0] );
f(:,:,4)=circshift( squeeze(f(:,:,4)), [+0,-1] );
f(:,:,5)=circshift( squeeze(f(:,:,5)), [+1,+1] );
f(:,:,6)=circshift( squeeze(f(:,:,6)), [-1,+1] );
f(:,:,7)=circshift( squeeze(f(:,:,7)), [-1,-1] );
f(:,:,8)=circshift( squeeze(f(:,:,8)), [+1,-1] );
end
% -------------------------------------------------------------------------------------
function f = bcs(M,N,u0,f)
% Left Wall(Bounce back)
f(1,:,1) = f(1,:,3);
f(1,:,5) = f(1,:,7);
f(1,:,8) = f(1,:,6);

% Bottom Wall(Bounce Back)
f(2:M,1,2) = f(2:M,1,4);
f(2:M,1,5) = f(2:M,1,7);
f(2:M,1,6) = f(2:M,1,8);

% Right Wall(Bounce Back)
f(M+1,:,3) = f(M+1,:,1);
f(M+1,:,6) = f(M+1,:,8);
f(M+1,:,7) = f(M+1,:,5);

rhon = zeros(M+1,N+1);
% Top Wall ( Moving with Lid Velocity u0)
for i = 2:M
    rhon(i,N+1) = f(i,N+1,9) + f(i,N+1,1) + f(i,N+1,3) +...
                    2*(f(i,N+1,2) + f(i,N+1,6) + f(i,N+1,5));
                
    f(i,N+1,4) = f(i,N+1,2);         
    f(i,N+1,7) = f(i,N+1,5) +0.5*(f(i,N+1,1)-f(i,N+1,3)) - (1/2)*rhon(i,N+1)*u0;
    f(i,N+1,8) = f(i,N+1,6) -0.5*(f(i,N+1,1)-f(i,N+1,3)) + (1/2)*rhon(i,N+1)*u0;
end
end

% -------------------------------------------------------------------------------
function[rho,u,v]=rhoU(M,N,f)
rho=sum (f,3);
for i=1:M+1
rho(i,N+1)=f(i,N+1,9)+f(i,N+1,1)+f(i,N+1,3)+2.*(f(i,N+1,2)+f(i,N+1,6)+f(i,N+1,5));
end
%calculate velocity compnents
u = ( sum(f(:,:,[1 5 8]),3) - sum(f(:,:,[3 6 7]),3) )./rho;
v = ( sum(f(:,:,[2 5 6]),3) - sum(f(:,:,[4 7 8]),3) )./rho;
end