The matter of the stability for multidimensional diffusion-advection-reaction problems treated with the semi-discretization method is remaining challenge because when all the stepsizes tend simultaneously to zero the involved size of the problem grows without bounds. Solution of such problems is constructed by starting with a semi-discretization approach followed by a full discretization using exponential time differencing and matrix quadrature rules. Analysis of the time variation of the numerical solution with respect to previous time level together with the use of logarithmic norm of matrices are the basis of the stability result. Sufficient stability conditions on stepsizes, that also guarantee positivity and boundedness of the solution, are found. Numerical examples in different fields prove its competitiveness with other relevant methods.