Reducing iteration computation time

Hi, I am trying to do finite difference calculation with MATLAB to estimate the total head within the flow region in a groundwater flow analysis. I have already done the code and it run well. Nonetheless, the computation time during the iteration procedure took a very long time for a very fine grid (nsize). Up to now, my guess is because of the size of the M matrix. The bigger the matrix is, the longer the calculation time. For nsize equals to 0.25, the dimension of M matrix is up to 97 x 217, and it required more than 30 mins to finish the computation of one model. If you can help with my problem, I really appreciate it. Thank you. The following code is the code for the iteration:

%%Step 3: Matrix iteration construction
M_mod = M;
while err > tol
  %Update iteration counter
  k = k+1;
  %Loop through computational nodes
  for i = 2:Ny;
    for j = 1:Nx; %Retained zone
      if j == 1;
        if i == Ny;
          M_mod(i,j) = 0.25*(2*M(i-1,j)+2*M(i,j+1));
        else
          M_mod(i,j) = 0.25*(M(i+1,j)+2*M(i,j+1)+M(i-1,j));
        end
      elseif j == (2*D/nsize); %Retained zone near sheeting pile
        if i <= ((He-dj+Hp)/nsize);
          M_mod(i,j) = 0.25*(M(i+1,j)+M(i-1,j)+2*M(i,j-1));
        elseif and(i > ((He-dj+Hp)/nsize),i <= Ny-1);
          M_mod(i,j) = 0.25*(M(i+1,j)+M(i-1,j)+M(i,j-1)+M(i,j+1));
        elseif i == Ny
          M_mod(i,j) = 0.25*(M(i,j+1)+2*M(i-1,j)+M(i,j-1));
        end
      elseif j == 2*D/nsize+1; %Nodes along the interface cut-off
        if and(i >= 2, i <=((He-dj+Hp)/nsize));
          M_mod(i,j) = M(i,j);
        elseif i == ((He-dj+Hp)/nsize+1);
          M_mod(i,j) = 0.25*(M(i,j+1)+2*M(i+1,j)+M(i,j-1));
        elseif i == Ny;
          M_mod(i,j) = 0.25*(M(i,j+1)+2*M(i-1,j)+M(i,j-1));
        else
          M_mod(i,j) = 0.25*(M(i+1,j)+M(i,j+1)+M(i-1,j)+M(i,j-1));
        end
      elseif and(and(j >= (2*D/nsize+2), j <= Nx),i <= ((He-dj+di)/nsize)+1); %Nodes in the excavation area
        M_mod(i,j) = M(i,j);
      elseif j == Nx; %Excavated zone
        if i == Ny;
          M_mod(i,j) = 0.25*(2*M(i-1,j)+2*M(i,j-1));
        else
          M_mod(i,j) = 0.25*(M(i+1,j)+2*M(i,j-1)+M(i-1,j));
        end
      elseif j == (2*D/nsize+2); %Excavated zone near sheeting pile
        if i == Ny;
          M_mod(i,j) = 0.25*(M(i,j+1)+2*M(i-1,j)+M(i,j-1));
        elseif and(i > ((He-dj+di)/nsize+1),i <= ((He-dj+Hp)/nsize));
          M_mod(i,j) = 0.25*(M(i+1,j)+M(i-1,j)+2*M(i,j+1));
        elseif and(i > ((He-dj+Hp)/nsize),i <= Ny-1);
          M_mod(i,j) = 0.25*(M(i+1,j)+M(i-1,j)+M(i,j+1)+M(i,j-1));
        end
     elseif i == Ny; %Global
       M_mod(i,j) = 0.25*(M(i,j+1)+2*M(i-1,j)+M(i,j-1));
     else
       M_mod(i,j) = 0.25*(M(i+1,j)+M(i,j+1)+M(i-1,j)+M(i,j-1));
     end
   end
 end
 %Calculate error
 err = sqrt(sum(sum((M_mod-M).^2)));
 %Update M_mod
 M = M_mod;
end