In chemical kinetics, the law of mass action gives us reaction rates of the form
$$r=k x^a y^b$$
where $r$ is the time derivative of either $x$ or $y$ times a constant
$$r=-\frac{dx}{\beta dt}=-\frac{dy}{\gamma dt}$$
Here's an example
$$r=k x^2 y= -\frac{dx}{2dt} = -\frac{dy}{dt}$$
In general, how would I go about solving these differential equations for $x(t)$ and $y(t)$?
On
If $x'=-\beta r$ and $y'=-\gamma r$ then $\gamma x'=\beta y'$ hence $\gamma x-\beta y$ is constant, say $$\beta y=\gamma (x-u_0)$$ with $$u_0=x_0-\beta \gamma^{-1}y_0$$ hence $$x'=-\beta r=-\beta k x^ay^b=-k_0x^a(x-u_0)^b$$ with $$k_0=\beta^{1-b}\gamma^bk$$ thus, $k_0>0$. The rest depends on the sign of $u_0$, thus, of the initial conditions $(x_0,y_0)$ and of the parameters $(\beta,\gamma)$, but basically, one translates this into the differential formulation $$\frac{dx}{x^a(x-u_0)^b}=-k_0dt$$ and then into its integrated version $$\int^{x_0}_{x(t)}\frac{du}{u^a(u-u_0)^b}=k_0t$$ which provides $x(t)$ as an implicit function of $t$ since, introducing any primitive $$g(u)=\int\frac{du}{u^a(u-u_0)^b}$$ defined in a neighborhood of $x_0$, one gets $$g(x(t))=g(x_0)-k_0t$$ that is,
$$x(t)=g^{-1}\left(g(x_0)-k_0t\right)$$
and the rest depends on how tractable the inverse function $g^{-1}$ is.
In the context of chemical kinetics, $a$ and $b$ are positive integers hence $$\frac1{u^a(u-u_0)^b}$$ is a rational fraction and, for every $c_0\ne0$, $g(u)$ is a linear combination of the functions $$1\qquad\log|u|\qquad\log|u-u_0|\qquad \frac1{u^k}\ (1\leqslant k\leqslant a-1)\qquad \frac1{(u-u_0)^\ell}\ (1\leqslant\ell\leqslant b-1)$$
Case $u_0=0$: This is when the initial conditions are such that $$\gamma x_0=\beta y_0$$ then one gets the primitive $$g(u)=\int\frac{du}{u^{a+b}}=-\frac1{cu^c}\qquad c=a+b-1\geqslant1$$ hence
$$x(t)=\frac{x_0}{(1+cx_0^ck_0t)^{1/c}}$$
Case $a=b=1$ and $u_0\ne0$: The solution is explicit as well since then, $$g(u)=-\frac1{u_0}\log\left|\frac{u}{u-u_0}\right|$$ hence $$\log\left|\frac{x(t)}{x(t)-u_0}\right|=\log\left|\frac{x_0}{x_0-u_0}\right|-u_0k_0t$$ that is, $$\frac{\gamma x(t)}{\gamma x(t)-\gamma x_0+\beta y_0}=\frac{\gamma x_0}{\beta y_0}\,\exp(-u_0k_0t)$$ which leads to
$$x(t)=\frac{x_0(\beta y_0-\gamma x_0)}{\beta y_0\exp((\beta y_0-\gamma x_0)kt)-\gamma x_0}$$
Numerically. In particular, if you have more than two reactions, pencil and paper methods will either fail or be too error prone.
If you do not want to do the programming yourself: standard software packages like octave or matlab have built-in solvers for systems of ordinary differential equations which are very accurate and fast.