For truly elementary reactions, the law of mass action can generally be assumed to hold fairly accurately, as predicted by the collisional theory of reaction kinetics. (This is something of a tautology, since a deviation from the law of mass action generally indicates that the reaction isn't really elementary.) In particular, this implies that:
- elementary unimolecular reactions occur at a rate proportional to the concentration of the reactant (or more generally to its activity, if we're treating e.g. the effects of a solvent on the reaction rate implicitly);
- elementary bimolecular reactions with distinct reactants occur at a rate proportional to the product of the concentrations (or activities) of the reactants;
- elementary reactions between two molecules of the same species occur at a rate proportional to the square of the concentration (or activity) of the reactant; and
- elementary reactions involving more than two reactants don't occur, except at negligible rates, since they would have to involve a perfectly simultaneous collision of three molecules. (Any such reactions that are observed to occur must actually proceed via one or more ephemeral transition states produced in bimolecular collisions, and are thus not actually elementary.)
Thus, in principle, any* such system of simultaneous elementary reactions can be modeled using a system of ordinary differential equations of the form:
$$
\frac{d}{dt} [A_x]
= \left( \sum_{y=1}^n k_{xy}[A_y] \right)
+ \left( \sum_{y=1}^n \sum_{z=1}^y k_{xyz}[A_y][A_z] \right),
$$
where $[A_1]$ to $[A_n]$ are the concentrations of the $n$ distinct chemical species present in the system (some of which might only exist as temporary transition states), $k_{xy}$ is the sum of the rate constants for reactions of the form $A_y \to A_x + (\text{other products})$, and $k_{xyz} = k_{xzy}$ is the sum of the rate constants for reactions of the form $A_y + A_z \to A_x + (\text{other products})$.**
So yes, if you know the rate constants for all the elementary mono- and bimolecular reactions your reactants can undergo, then you can indeed plug them into a differential equation like the one above and solve it.
However, as a fair warning from someone who's worked on numerically solving these kinds of differential equation systems, they are notorious for often exhibiting high stiffness and consequent numerical instability. This is particularly often the case if some of the reaction constants are much larger and/or some of the reactant concentrations are much higher than others. Efficiently solving such stiff ODE systems is a hard problem, and there's a whole bunch of advanced (and often jargon-filled) mathematical literature on it.
In the end, the practical ways of solving such systems usually come down to timescale separation or other stiffness-reducing approximations, where we essentially assume that certain "fast" reactions occur so rapidly that the concentrations of their products and reactants are always at a local quasi-equilibrium, where the overall net reaction rate of the fast reactions and their reverse reactions is zero. This allows us to then only focus on the remaining "slow" reactions, which are presumably the ones we're actually interested in.
Alas, as a side effect, these approximations also often tend to destroy the mathematical simplicity of the basic elementary reaction model by reintroducing non-elementary reaction laws and/or nontrivial dependencies between the reactant concentrations. Also, unfortunately, at least I know of no general method for automatically constructing such approximations, just as there's no fully general method for efficiently solving the original stiff ODE systems in the first place.
*) At least assuming that the reactions occur in a homogeneous well-mixed environment. Things like reactions occurring on the surface of a substrate can be included by treating the substrate-bound and free forms of a molecule as distinct species and, if necessary, also treating free surface bonding sites as a separate species with its own concentration. Similar tricks can also be used to model e.g. reactions occurring in two immiscible solvents, provided that mixing within each solvent phase can be assumed to be rapid. However, if the reactant concentrations may vary spatially to a significant degree, then one would have to use a reaction–diffusion model instead.
**) If the products of a reaction include $m$ molecules of species $A_x$, that reaction's rate constant should be included $m$ times in $k_{xy}$ or $k_{xyz}$. Also note that, to avoid double counting $A_y + A_z \to A_x + \dots$ and $A_z + A_y \to A_x + \dots$ as distinct reactions, the last sum in the differential equation is only taken over $z \le y$.
