[SashaPrival]

Hello. Newbie here. I am following this guide to create my very first simple CFD solvers for incompressible 2D fluid in a lid-driven square cavity. In the guide, the pressure is considered with zero-derivative Neumann BCs imposed on all boundaries. I have managed to write and run a Matlab code following the guide. It runs fine.

Then I wanted to change the top boundary condition for the pressure to 0 Dirichlet like in the Step 11 of famous beginner guide of Lorena Barba. I have modified the Laplacian matrix contents to include top Dirchlet BC. Assigned 0 to the corresponding elements of the source (right-hand side) vector. Forced 0 value to the corresponding elements of the solution (pressure) vector. Here is the complete Matlab code I wrote:

clear; clc;
% grid parameters
Lx=1; % x length
Ly=1; % y length
nx=41; % number of grid points along x
ny=41; % number of grid points along y

%simulation parameters
dt=0.001; % time step
t_final=5.0; % simulation run time

% Number of timesteps
Nt=t_final/dt;

% physics parameters
rho=1.0; % density 
nu=0.001; % viscosity coefficient (kinematic)

% Index extents
imin=2;
imax=imin+nx-1; % ok
jmin=2;         % ok
jmax=jmin+ny-1;

%create mesh sizes
dx=Lx/(nx);
dy=Ly/(ny);
dxi=1/dx;
dyi=1/dy;

% preallocate important arrays 
p=zeros(imax,jmax); % pressure array

%p=zeros(nx,ny); % ok
us=zeros(imax+1,jmax+1); % u - velocity without pressure
vs=zeros(imax+1,jmax+1); % v - velocity without pressure
R = zeros(nx*ny, 1);  % RHS of Pressure Poisson Equation
u=zeros(imax+1,jmax+1); % corrected u
v=zeros(imax+1,jmax+1); % corrected v
L=zeros(nx*ny,nx*ny); % Laplacian matrix

% initial values
t=0; % initial time
u_bot=0; % initial velocity for bottom wall
u_top=1; % initial velocity for top wall
v_lef=0; % initial velocity for left wall
v_rig=0; % initial velocity for right wall

% Creat Laplacian operator for solving pressure Poisson equation
for j=1:ny %number of block matrices along y
    for i=1:nx %number of block matrices along x
        q = i+(j-1)*nx;
        %if(mod(q,1*nx)==0)
        if (q > (nx-1)*ny) % if top boundary            
            L(q,q)=1; % Zero Dirichlet condition
        else
            L(q,q)=-(2*dxi^2+2*dyi^2);
            for ii=i-1:2:i+1
                if (ii>0 && ii<=nx) % Interior points
                    L(q,ii+(j-1)*nx)=dxi^2;      
                else % Neuman conditions on boundary
                    L(q,q)= L(q,q)+dxi^2;
                end
            end
            for jj=j-1:2:j+1
                if (jj>0 && jj<=ny) % Interior points
                    L(q,i+(jj-1)*nx)=dyi^2;
                else % Neuman conditions on boundary
                    L(q,q)= L(q,q)+dyi^2;
                end
            end
        end
    end
end
%disp("L="); disp(num2str(L));

% solver
disp("Starting simulation ****************************");


while t <= t_final
    % update time
    t = t + dt;
    clc; disp("t:"); disp(t);


    u_top=1;


    % u velocity predictor
    for j = jmin:jmax
        for i = imin+1:imax
            v_here = 0.25*(v(i-1,j)+v(i-1,j+1)+v(i,j)+v(i,j+1));
            us(i,j)=u(i,j)+dt* ...
                (nu*(u(i-1,j)-2*u(i,j)+u(i+1,j))*dxi^2 ...
                +nu*(u(i,j-1)-2*u(i,j)+u(i,j+1))*dyi^2 ...
                -u(i,j)*(u(i+1,j)-u(i-1,j))*0.5*dxi ...
                -v_here*(u(i,j+1)-u(i,j-1))*0.5*dyi);
        end
     end


    % v velocity predictor
    for j = jmin+1:jmax
        for i = imin:imax
            u_here = 0.25*(u(i,j-1)+u(i+1,j-1)+u(i,j)+u(i+1,j));
            vs(i,j)=v(i,j)+dt* ...
                (nu*(v(i-1,j)-2*v(i,j)+v(i+1,j))*dxi^2 ...
                +nu*(v(i,j-1)-2*v(i,j)+v(i,j+1))*dyi^2 ...
                -u_here*(v(i+1,j)-v(i-1,j))*0.5*dxi...
                -v(i,j)*(v(i,j+1)-v(i,j-1))*0.5*dyi);
        end
    end


    % compute R.H.S. (R) of the PPEquation using vs & us.
    n=0;    
    %R = zeros(nx*ny, 1);  % RHS of PPEquation
    for j=jmin:jmax
        for i=imin:imax
            n=n+1;
            R(n)=(rho/dt)*((us(i+1,j)-us(i,j))*dxi+(vs(i,j+1)-vs(i,j))*dyi);
            % if(mod(n,nx)==0)
            if (n > (nx-1)*ny) % if top boundary
                R(n)=0; % force 0 value
            end
            %disp(["R[",n,"]:",R(n)]);
        end
     end


    % Find pressure solution vector using direct method
    pv=L\R;
    % disp("pv="); disp(pv);


    n=0;    
    p=zeros(imax,jmax);
    % convert pressure vector into mesh array
    for j=jmin:jmax
       for i=imin:imax
           n=n+1;
           % if (mod(n,nx)==0)
           if(n > (nx-1)*ny) % if top boundary
               pv(n)=0; % force Zero Dirichlet
           end
           p(i,j)=pv(n);
       end
    end
    
    %disp("p="); disp(num2str(p));

    % correcting (updating) u & v velocities


    for j=jmin:jmax
        for i=imin+1:imax
            u(i,j)=us(i,j)-(dt/rho)*(p(i,j)-p(i-1,j))*dxi;
        end
    end


    for j=jmin+1:jmax
        for i=imin:imax
            v(i,j)=vs(i,j)-(dt/rho)*(p(i,j)-p(i,j-1))*dyi;
        end
    end

    % set Boundary Conditions
    u(:,jmin-1)=u(:,jmin)-2*(u(:,jmin)-u_bot);
    u(:,jmax+1)=u(:,jmax)-2*(u(:,jmax)-u_top);
    v(imin-1,:)=v(imin,:)-2*(v(imin,:)-v_lef);
    v(imax+1,:)=v(imax,:)-2*(v(imax,:)-v_rig);
    
end;


disp("End of simulation");
disp("Last t="); disp(num2str(t));
disp("Last p="); disp(num2str(p));
disp("Last u="); disp(num2str(u));
disp("Last v="); disp(num2str(v));
disp("Finita!");

But it leads to NAN after some time of running (~3.1 s). I tried to thoroughly scan my code for bugs or other types of errors many times. But I can't find anything. I tried to decrease time step (from 1e-2 to 1e-4). Tried to increase grid resolution (from 11x11 to 51x51). No luck. Maybe there is some fundamental flaw in my realization?
P.S. This is my first post here. Please kindly remind me of any "don'ts" I have made in the post if there are any. Thank you.

[SashaPrival]

Is there someone who will give me some advice or help in anyway?

[zacko]

Did you check at which loop it crashes? Did you print out your matrices? Maybe check the initialalization first and iterate once and recheck the flow field. Maybe there are any obvious problems already apparent