Integrals on $[0,\infty)$ where the integrand is of the form $Q^m(a\sqrt{x})\,{\rm p}(x)$, where $Q$ is the Gaussian $Q$ function, ${\rm p}(\cdot)$ a Gamma PDF, $m$ a positive integer and $a>0$; or of the form $\erf^m(ax+b)\,x^n\exp(-c^2x^2+2dx)$, where $\erf(x)$ is the error function, and $n$ a non-negative integer, arise in performance modelling of communication and machine learning systems. Such integrals cannot be evaluated analytically in general, but they are reducible to a key integral whose integrand is $\erf(ax+b)\,{\rm N}(x;m,s)$ where ${\rm N}()$ is a Gaussian PDF with mean $m$ and variance $s$. Seeking an efficient and accurate evaluation method, we develop a new 4-term exponential quadratic approximator (EQA) for the error function that includes both linear and quadratic terms in its exponents. The EQA minimises a sum-of-squares cost function with two ``spline-type'' constraints, i.e., constraints on the function value and its first derivative. This 12-parameter constrained optimisation problem is reduced to an unconstrained 8-parameter one by inverting a 4-D linear system, then solved by gradient descent. The resulting approximator has a maximum absolute error of $1.65\times 10^{-4}$ on the real line, and outperforms many other exponential sum approximators for $\erf(x)$ on $x\in[0,1.5]$ and for $Q(x)$ on $x\in[0,2]$. Moreover, due to its functional form, the EQA leads to an analytical formula for the key integral that is accurate to 3 to 4 significant figures when compared with known special cases or computationally intensive Monte Carlo integration. The EQA can equally be used to obtain closed forms for the average symbol error probability of various modulation schemes on Rayleigh fading channels.