Bifurcation theory is a powerful tool for studying how the dynamics of a neural network model depends on its underlying neurophysiological parameters. However, bifurcation theory of neural networks has been developed mostly for mean-field limits of infinite-size spin-glass models, for finite-size dynamical systems whose units have a graded, continuous output, and for models with discrete-output neurons that evolve in continuous time. To allow progress on understanding the dynamics of some widely used classes of neural network models with discrete units and finite size, which could not be studied thoroughly with the previous methodology, here we introduced algorithms that perform a semianalytical bifurcation analysis of a finite-size firing-rate neural network model with binary firing rates and discrete-time evolution. In particular, we focus on the case of small networks composed of tens of neurons, to which existing statistical methods are not applicable. Our approach is based on a numerical brute-force search of the stationary and oscillatory solutions of the model, from which we derive analytical expressions of its bifurcation structure by means of the state-to-state transition probability matrix. Our algorithms determine how the network parameters affect the degree of multistability, the emergence and the period of the neural oscillations, and the formation of spontaneous symmetry breaking in the neural populations. While this technique can be applied to networks with arbitrary (generally asymmetric) connectivity matrices, in particular we introduce a highly efficient algorithm for the bifurcation analysis of sparse networks. We also provide some examples of the obtained bifurcation diagrams and a python implementation of the algorithms.