first-order binding kinetics with multiple identical binding sites
The scheme you need seems to be of the form 
$$\ce{S + M_i  &lt;=&gt;[k_2][k_{-1}] M_{i+1}    }$$
where S is the molecule in solution and M the bead with i sites filled. 
If there is no restriction on how many S may be associated with any bead and as the  S are indistinguishable molecules the number of ways of arranging the S among a total of M beads is 
$$ W=\frac{M!}{\prod M_i!}$$
and the constrains are $\Sigma_0^N M_i = M$ and $\Sigma iM_i = N$ where N is the total number of S molecules.
Solving with these constraints is done using Lagrange multipliers, in the same manner as the Statistical Mechanics calculation of Boltzmann distribution, and leads to the distribution of finding i molecules of type S in a bead as,
$$p_i = \frac{\bar n ^i}{(1+\bar n)^{i+1}}$$
where $\bar n$ is the average occupancy of a bead. This distribution, which is clearly not Poissonian, has the form that it decreases as i increases even for moderate $\bar n$.

If the kinetic scheme is changed slightly to 
$$\ce{S + M_i  &lt;=&gt;[k_2][(i + 1)k_{-1}] M_{i+1}    }$$
where the ‘off’ rate constant now depends on how many S are on a bead then the probability becomes Poissonian. This change/trick effectively makes the molecules distinguishable and now
$$p_i = \frac{\bar n^i}{i!}e^{-\bar n}$$
with $\bar n = k_2\ce{[S]}/k_{-1}$.

Which of these or other models applies will have to be determined by experiment.

If you want to solve the equations numerically then it is best to try a Master equation approach.

The set of equations for the population of M is written as 
$$d\vec M/dt = \vec k \vec M$$
and has the formal solution $\vec M=\vec M_0e^{\vec kt}$ where $\vec M$ is a vector of concentrations of species i and $\vec k$ a matrix of rate constants such as

$$ \bar k=\left [\begin {matrix} 
-Sk_2 &amp; 2k_{-1} &amp; 0 &amp; 0 &amp; 0 &amp; \cdots\\
Sk_2 &amp; -Sk_2 -2k_{-1} &amp; 3k_{-1} &amp; 0 &amp; 0 &amp; \cdots \\
\vdots&amp; \vdots&amp;\vdots &amp; \vdots&amp; \vdots&amp; \\
\end {matrix}\, \right ]$$

where S is the concentration of S and is assumed to be in excess so that its value is constant. (If you want S to vary then you will need to do a numerical integration on all N equations instead of this method).

To solve, the secular determinant is set up and the eigenvalues $\lambda$ obtained. The populations at time t are given by
$$M(t)= x[e^{\lambda t}]x^{-1}M_0$$
where x is the matrix of eigenvectors corresponding to eigenvalues $\lambda$ and $[e^{\lambda t}] $ is a diagonal matrix of $e^{\lambda_1t}, e^{\lambda_2 t} \dots $ etc.