Numerical methods for solving kinetic equations for diffusion-influenced bimolecular reactions are presented for three cases. Finite difference method is used to solve diffusion-reaction equations for the pair distribution function. The kinetic equation for the concentration is evolved by the Runge-Kutta method with adaptive time step. The boundary doubling method is introduced to study long time dynamics, where the truncation problem of the infinite boundary is crucial. The above methods are applied, in the first case, to the classical Smoluchowski approach to a binary reaction with random initial condition and the results are compared with ones in two dimension. In the second case, an isolated pair recombination dynamics with a delta function initial condition is investigated and the results are compared with analytic expression in three dimension with spherical symmetry. A more complicated system with the hierarchical Smoluchowski approach with the Kirkwood superposition approximation is also investigated in the third case. The efficiency and the accuracy of the numerical calculations are examined against the asymptotic analytical solutions and a Monte Carlo simulation in one dimension.