MacKay Neal Algorithm - Stuffing a matrix with ones and zeros

I need to fill in a matrix with ones and zeros. The parameters that are provided are its dimension (N columns and k = N * (1-r) rows) and two vectors h and v that represent the distribution of rows and columns, respectively. For example, suppose we have a matrix of dimension 5x10, whose vectors h and v take the following values:

v = [0 0.2 0 0.4 0.4];
h = [0 0.1 0.2 0 0.1 0.3 0 0 0.2 0.1];

This means that we have the following information about the rows of H:

• 10% of the total contains 2 ones (and the rest of their elements are zeros).
• 10% of the total contains 5 ones (and the rest of their elements are zeros).
• 10% of the total contains 10 ones.
• 20% of the total contains 3 ones (and the rest of their elements are zeros).
• 20% of the total contains 9 ones (and the rest of its elements -in this case, one- are zeros).
• 30% of the total contains 6 ones (and the rest of their elements are zeros).
• None of its rows has a single, or 4, or 7, or 8, (or none) one.

And, as for the H columns:

• 20% of the total of them contains 2 ones (and the rest of their elements are zeros).
• 40% of the total of them contains 4 ones (and the rest of their elements -in this case, one- are zeros).
• 40% of the total of them contains 5 ones.
• None of its columns has a single, or 3, (or none) one.

Under these premises, I used the following code:

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% MacKayNeal Function                                                     %
%                                                                         %
% Inputs                                                                  %
% N: length of the codeword (Number of columns of parity-check matrix H)  %
% r: rate of the LDPC code                                                %
% v: distribution of edges                                                %
% h: distribution of nodes                                                %
%                                                                         %
% Outputs                                                                 %
% H: parity-check matrix                                                  %
%                                                                         %
% Author: Aitor López Hernández                                           %
%                                                                         %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

function H = MacKayNeal(N,r,v,h)

H = zeros(N*(1-r),N); % Creation of an idle parity-check matrix

a = []; % Number of 1s per column of H
for i=1:length(v)
    for j=1:(v(i)*N)
        a = [a,i];
    end
end

b = []; % Number of 1s per row of H
for i=1:length(h)
    for j=1:(h(i)*(N*(1-r)))
        b = [b,i];
    end
end

b_i = 1:length(b); % Indexes of rows of H that still have 1s to append

for i=1:N
    for j=1:1:length(b)
    if b(j) == 0
        b_i = setdiff(b_i,b(j));
    end
    end
    c = datasample(b_i,a(i),'Replace',false);
    c_i = zeros(1,length(b));
    for j=1:a(i) 
        H(c(j),i)=1;
        c_i(1,c(j))=1;
    end
    b = b - c_i;
end

% Removal of length-4 cycles

for i=1:N*(1-r)-1
    for j=i+1:N*(1-r)
        w = and(H(i,:),H(j,:));
        c1 = find(w);
        lc = length(c1);
        if lc > 1
            % If found, flip one 1 to 0 in the row with less number of 1s
            if length(find(H(i,:))) < length(find(H(j,:)))
            % Repeat the process until only one column left
            for cc = 1:(lc-1)
                 H(j,c1(cc)) = 0;
            end
            else
            for cc = 1:(lc-1)
                H(i,c1(cc)) = 0;
            end
            end
        end
    end
end
end

(The part following the comment "Removal of 4-length cycles" is not relevant to my query).

The problem lies in the multiplications that appear in the loops in lines 5 and 11. Their results must necessarily be integers for the algorithm to work (keeping the sum of the elements of v and h equal to 1, since its elements represent probability distributions).

The following values can be used as an example to show the malfunction of my code. Although the program runs without errors, it is observed that it does not do it properly, since the vector b does not have a size of N (1-r) = 50 but of 48 elements.

N = 100;
r = 0.5;
h = [0 0.17 0.03 0.2 0 0 0.17 0.3 0 0 0 0 0.03 0.1];
v = [0 0.1 0 0 0 0.7 0.2];

This algorithm needs to work for values much greater than N = 100 and other factors other than ar = 0.5 -and the same applies to a and h, which could also have values with many more decimals (since they are obtained as a result of a problem of optimization with the 'fmincon' function).

Does anyone have a solution for my problem?

Thank you very much in advance, and greetings.