The present study focuses on the numerical solution of the transient convection-diffusion-reaction equation by transforming it into modified Helmholtz equation through an exponential type transformation. In the spatial discretization of the problem domain two different boundary element methods (BEM), namely the domain BEM (DBEM) and the dual reciprocity BEM (DRBEM), are employed which are combined with an implicit backward finite difference time integration. The BEM techniques differ in the sense of treating the leftover domain integral. That is, in DBEM the domain integral is kept and calculated by quadrature while it is reduced to an equivalent boundary interal by means of radial basis functions in DRBEM. The numerical simulations are first carried out for several values of diffusion coefficient in the transient convection-diffusion-reaction equation. The results reveal that DBEM gives more accurate results for smaller values of diffusion coefficient compared to DRBEM. Thus, DBEM is further used for the solution of the coupled transient convection-diffusion type magnetohydrodynamic flow equations. The results are presented by equivelocity and current lines for various values of problem parameters, which show the well-known characteristics of the MHD flow.