Optimize algorithm that generates the number of units in each binary state

Question

TL;DR: I need to find all possible combinations of N row vectors (of size 1xB), whose row-wise sum produces the desired result vector (also of size 1xB).

---

I have a binary matrix (1 or 0 entries only) of size N x B where N denotes the number of units and B denotes the number of bins. Each unit, i.e., each row, of the matrix can be in one of 2^B states. That is, if B=2, the states possible are {0,0}, {0,1}, {1,0} or {1,1}. If B=3, then the possible states are {0,0,0}, {0,0,1}, {0,1,0}, {0,1,1}, {1,0,0}, {1,0,1}, {1,1,0} or {1,1,1}. Basically the binary representation of the numbers from 0 to 2^B-1.

For the matrix, I know the sum over the rows of the matrix, for example, {1,2}. This sum can be achieved through different binary matrices like [0,0;0,1;1,1] or [0,1;0,1;1,0]. The number of units in each state are {1,1,0,1} and {0,2,1,0}, respectively for each of the matrices, where the first number corresponds to the first state {0,0}, second to the second state {0,1} and so on in increasing order. My problem is to find all possible vectors of these numbers of states that satisfy a particular matrix sum.

Now to implement this in MATLAB, I used recursion and a global variable. This to me was the easiest approach, however, it takes a lot of time. The code I used is given below:

function output = getallstate()
   global nState % stores all the possible vectors
   global nStateRow % stores the current row of the vector
   global statebin %stores the binary representation of all the possible states
   nState = [];
   nStateRow = 1;
   nBin = 2; % number of columns or B
   v = [1 2]; % should always be of the size 1 x nBin
   N = 3; % number of units
   statebin = de2bi(0:(2 ^ nBin - 1), nBin) == 1; % stored as logical because I use it to index later
   getnstate(v, 2 ^ nBin - 1, nBin) % the main function
   checkresult(v, nState, nBin) % will result in false if even one of the results is incorrect
   % adjust for max number of units, because the total of each row cannot exceed this number.
   output = nState(1:end-1, :); % last row is always repeated (needs to be fixed somehow)
   output(:, 1) = N - sum(output(:, 2:end), 2); % the first column, that is the number of units in the all 0 state is always determined by the number of units in the other states
   if any(output(:, 1) < 0)
      output(output(:, 1) < 0, :) = [];
   end
end

function getnstate(r, state, nBin)
   global nState
   global nStateRow
   global statebin

   if state == 0
      if all(r == 0)
         nStateRow = nStateRow + 1;
         nState(nStateRow, :) = nState(nStateRow - 1, :);
      end
   else
      for a = 0:min(r(statebin(state + 1, :)))
         nState(nStateRow, state + 1) = a;
         getnstate(r - a * statebin(state + 1, :), state - 1, nBin);
      end
   end
end

function allOk = checkresult(r, nState, nBin)
   % just a function that checks whether the obtained vectors all result in the correct sum
   allstate = de2bi(0:(2 ^ nBin - 1), nBin);
   allOk = true;
   for iRow = 1:size(nState, 1)
      sumR = sum(bsxfun(@times, allstate, nState(iRow, :).'), 1);
      allOk = allOk & isequal(sumR,r);
   end
end

function b = de2bi(d, n)
d = d(:);
[~, e] = log2(max(d));
b = rem(floor(d * pow2(1-max(n, e):0)), 2);
end

The above code works fine and gives all possible states but, as is expected, it gets slower as you increase the number of columns (B) and the number of units (N). Also, it uses globals. The following are my questions:

1. Is there a way to generate these without using globals?
2. Is there a non-recursive way for this algorithm?

EDIT 1

3. In what way do the above and still have an optimised algorithm which is faster than the current version?

EDIT 2

Added the de2bi function to remove dependency on the Communications Toolbox.