Working in the many-boson formalism is not that different than the more standard "gate-based" approach you might be more used to. The main difference is that the unitary operations technically act on infinite-dimensional spaces, because each mode can accomodate an arbitrary number of bosons.
However, as long as you stick with linear operations, meaning in this context operations that preserve the total number of bosons, things can be again represented as just matrices applied to vectors. So, for example, the $U_i$ operations you seem to be defining here can be understood as just standard beamsplitter operations with elements
$$\begin{pmatrix}\sqrt{1-R} & i\sqrt R \\ i\sqrt R & \sqrt{1-R}\end{pmatrix},$$
wrt the appropriate basis. Of course, the major catch is that this is only true when these operations are applied to input single bosons. The moment you inject, say, one boson in each input mode of the operation, this naive matrix description fails, and you get e.g. HOM-like effects.
You can still describe each operation in a linear optical context as a matrix in the above sense, even when more bosons are used, but you'll have to be more careful about the associated matrix representation. For example, the action of an $n$-mode unitary $U$ acting on $k$-boson input states will correspond to a matrix $\varphi_k(U)$ of dimension $\binom{n+k-1}{k}$, whose elements correspond, modulo factors, to permanents of $k\times k$ submatrices of the "standard" matrix representation of $U$. Though explaining how this works precisely might be a bit OT here.
Alternatively, you can of course just use the formalism with creation/annihilation operators to carry out the calculations. To find the "post-selected states" you then just need to find the coefficients of the terms containing the correct numbers of creation operators in the modes you want them to be in.
