Small relative perturbations to the entries of an essentially nonnegative matrix introduce small relative errors to entries of its exponential. It is thus desirable to compute the exponential with high componentwise relative accuracy. Taylor series approximation coupled with scaling and squaring is used to compute the exponential of an essentially nonnegative matrix. An a priori componentwise relative error bound of truncation is established, from which one can choose the degree of Taylor series expansion and the scale factor so that the exponential is computed with desired componentwise relative accuracy. To reduce the computational cost, the degree of the Taylor series expansion is chosen small, while the scale factor is chosen sufficiently large to achieve the desired accuracy. The rounding errors in the squaring stage are not serious as squaring is forward stable for nonnegative matrices. We also establish a posteriori componentwise error bounds and derive a novel interval algorithm for the matrix exponential of an essentially nonnegative matrix. Rounding error analysis and numerical experiments demonstrate the efficiency and accuracy of the proposed methods.