I had presented the program in pieces. The full program is below with the pieces in place.

That program assumes the queens are all in separate quadrants of the board. A variation of a section:

for q1r=1:4, for q1c=1:4

  queen(1,1)=q1r; queen(1,2)=q1c;

for q2r=1:4, for q2c=1:4

  queen(2,1)=q2r; queen(2,2)=q2c;

for q3r=5:8, for q3c=1:4

  queen(3,1)=q3r; queen(3,2)=q3c;

for q4r=1:8, for q4c=5:8

  queen(4,1)=q4r; queen(4,2)=q4c;

  

places two queens in one quadrant while the other two must still be in separate quadrants. It covers all such cases, but doesn't attempt to find all rotations and reflections, but that doesn't matter as such things are filtered out anyway.

This variation doesn't find any solutions, and the fewest locations unchecked by the queens was 3, rather than 2.

clearvars,clc

goodGrids={};

bestCt=999;

bestGrid=zeros(8);

for q1r=1:4, for q1c=1:4

  queen(1,1)=q1r; queen(1,2)=q1c;

for q2r=1:4, for q2c=5:8

  queen(2,1)=q2r; queen(2,2)=q2c;

for q3r=5:8, for q3c=1:4

  queen(3,1)=q3r; queen(3,2)=q3c;

for q4r=5:8, for q4c=5:8

  queen(4,1)=q4r; queen(4,2)=q4c;

  

  grid=zeros(8); g0=grid;

  

  for i=1:4

    g0(queen(i,1),queen(i,2))=1;

    diff=queen(i,2)-queen(i,1); diff2 = 9-queen(i,1)-queen(i,2);

    grid=grid+(diag(repmat(1,length(grid)-abs(diff),1),diff));

    grid=grid+flip(diag(repmat(1,length(grid)-abs(diff2),1),-diff2)); 

    grid(queen(i,1),:)=grid(queen(i,1),:)+1;

    grid(:,queen(i,2))=grid(:,queen(i,2))+1;

  end

  gct=grid==0;

  ct0=sum(gct(:));

  if ct0<bestCt

    bestCt=ct0;

    bestGrid=grid;

    goodGrids={grid g0};

  elseif ct0==bestCt

    good=true; gridA=grid;

    for i=1:5

      for j=1:length(goodGrids)

        if isequal(gridA,goodGrids{j})

          good=false;

          break

        end        

      end

      gridA=rot90(gridA);

    end

    gridA=flip(gridA);

    for i=1:5

      for j=1:length(goodGrids)

        if isequal(gridA,goodGrids{j})

          good=false;

          break

        end        

      end

      gridA=rot90(gridA);

    end

    if good

      goodGrids=[goodGrids  {grid g0}];

    end

  end

end

end

end    

end

end

end

end    

end

bestCt

length(goodGrids)/2

% goodGrids{:}

for i=1:2:length(goodGrids)

  grid1=goodGrids{i}; grid2=goodGrids{i+1};

  grid3=repmat("-",8);

  grid3(find(grid2==1))="Q";

  grid3(find(grid1==0))="b"

end

function [r,c]=rowcol(lin)

  r=floor(lin/8-1)+1;

  c=mod(lin/8-1,8)+1;

end