Solving the simplest coupled nonlinear ODES for chemical kinetics I am just trying to get the integrated form for the kinetics of the reaction 
$A + B \rightarrow C + D$ characterized by: 
$$
 -\dfrac{d[A]}{dt} = -\dfrac{d[B]}{dt} = k[A][B]  \;  .
$$
As you note, this is a system of two nonlinear coupled ordinary differential equations, I will put that in the familiar notation for mathematicians or physicists:
$$
\dfrac{dx(t)}{dt} = -k x(t) y(t)
$$
$$
\dfrac{dy(t)}{dt} = -k x(t) y(t)
$$
Can anyone helping me with some reference or ideas to solve it? 
 A: Since the rate of change of $ x $ is the same as the rate of change of $y $ you really only a single equation of with one variable. We write,
\begin{equation} 
x = y + c
\end{equation} 
where the constant $ c $ is determined by your initial conditions,
\begin{equation} 
c = x  (0) - y (0) 
\end{equation} 
(in your case it is the difference between the concentrations at the start). Inserting in this relation we have,
\begin{equation} 
x' = -k (  x ^2 - x c  ) 
\end{equation} 
This is now a simple equation to solve,
\begin{equation} 
\int _{x _0 } ^{ x} \frac{ d x }{ - k ( x ^2 - x c ) } = \int _0 ^t d t 
\end{equation} 
where $ x _0 \equiv x ( 0 ) $. Resorting to Mathematica (though you could use partial fractions) gives,
\begin{equation} 
   \log \left[ \frac{ x }{ x _0 } \frac{ x _0 - c }{ x - c } \right]    = kt
\end{equation} 
Isolating for $ x $ I get,
\begin{equation} 
x = - \frac{ c e ^{ kt} }{ \alpha - e ^{ kt} }
\end{equation} 
where $ \alpha \equiv ( x _0 - c ) / x _0 $. 
Plotting $x$[red] and $y$[blue] for different $x0$ and setting $c=1$ we have,
$\hspace{2cm}$
A: Hints: 


*

*Conclude that $y-x=c$ is a constant.

*Use separation of variables $-k\int \!\mathrm{d}t= \int \!\frac{\mathrm{d}x}{x(x+c)}$.
A: The other answers address how to solve this analytically, but I like numerical solutions to things so here goes:
$$\frac{d}{dt} \begin{bmatrix}x \\ y \end{bmatrix} = -\begin{bmatrix}kxy\\kxy\end{bmatrix}$$
which can be solved using any number of numerical methods. For simplicity, we can take the second-order Runge-Kutta method where $i$ is the time index. 
Step 1:
$$\begin{aligned}x^{i+1/2} &= x^i - \Delta t k x^i y^i\\y^{i+1/2} &= y^i - \Delta t k x^i y^i\end{aligned}$$
Step 2:
$$\begin{aligned}x^{i+1} &= \frac{1}{2}\left(x^i + x^{i+1/2} - \Delta t k x^{i+1/2} y^{i+1/2}\right)\\y^{i+1} &= \frac{1}{2}\left(y^i + y^{i+1/2}- \Delta t k x^{i+1/2} y^{i+1/2}\right)\end{aligned}$$
This method is then marched until you reach your physical time of interest or until you reach steady state (which is $t \rightarrow \infty$ technically but may be truncated by measuring the residuals and stopping when they approach zero.
Some care must be taken when choosing your time step, $\Delta t$. Too large and the integration will be unstable. Too small and you're wasting time. 
Analytical methods are awesome and all but the numerical approach will work for much more complicated mechanisms and with way more coupled equations. So it's a good tool to have in the box.
