I have been in many such instances and one such instance arose in consideration of the wavepacket size effects on neutron time-of-flight spectroscopy, circa 2009. Calculation of the expected time-of-flight of a Gaussian wave-packet led to an integral which is a specific case of the integral,
\begin{equation}\label{integral}
\int_0^{\infty} e^{-a^2 x^2} \mathrm{erf}(b x) \sin(s x)\, dx ,
\end{equation} where $\mathrm{erf}(z)$ is the error function, and $a$, $b$ and $s$ are complex parameters. GR, Mathematica and Maple failed me big time on this integral at the time. And so I had to evaluate it my way and the result is
\begin{equation}\label{sinintegral}
\int_0^{\infty} e^{-a^2 x^2} \mathrm{erf}(b x) \sin(s x)\, dx = \frac{\sqrt{\pi}}{2 a} e^{-s^2/4 a^2} \mathrm{erfi}\!\left(\frac{b s}{2 a \sqrt{a^2+b^2}}\right),
\end{equation} where $\mathrm{erfi}(z)$ is the imaginary error function and the equality holds for all complex $s$, $a$ and $b$ with $a\neq 0$ and $\mathrm{Re}(a^2+b^2)>0$.
More than eight years later, the latest iterations of Mathematica (version 11.3) and Maple (version 2018.2) still return the integral unevaluated even for specific values of the parameters, and GR’s $8^{\mathrm{th}}$ edition (edited by D. Zwillinger and V. Moll) still does not tabulate the integral. So let us establish the above integral and in the process learn some new tricks that we can keep in our bag of math hacks.
The primary trick is to show that the integral satisfies a differential equation and that its value is a particular solution to the differential equation determined by an initial condition. The differential equation is discovered by means of the classic method of integration by parts. To proceed, let us, in the meantime, limit ourselves to some fixed non-zero real $a$ and $b$, and to a real variable $s$, and define the function
\begin{equation}\label{function}
f(s)=\int_0^{\infty} e^{-a^2 x^2} \mathrm{erf}(b x) \sin(s x)\, dx .
\end{equation} The integral exist for all real $s$ so that $f(s)$ takes the entire real line as its domain; moreover, $f(s)$ is infinitely differentiable. The function $f(s)$ is exactly the integral \eqref{integral} only that its parameters and variable are real valued. Later we will appeal to analytic continuation to extend our result to complex $a$, $b$ and $s$. Now we perform integration by parts on the right hand side of equation \eqref{function} with $dv=\sin(s x) dx$ and $u=e^{-a^2 x^2} \mathrm{erf}(b x)$ to yield
\begin{eqnarray}
&&f(s)=\frac{2 b}{\sqrt{\pi} s}\int_0^{\infty} e^{-(a^2+b^2) x^2} \cos(s x)\, dx \nonumber\\
&&\hspace{20mm}-\frac{2 a^2}{s} \int_0^{\infty} e^{-a^2 x^2} x \, \mathrm{erf}(bx) \cos(s x) dx ,\label{byparts}
\end{eqnarray} where we have used the fact that the derivative of $\mathrm{erf}(x)$ is equal to $2 e^{-x^2}/\sqrt{\pi}$.
Under the given condition, the integral in the first term of equation \eqref{byparts} exists and is tabulated, and is given by\begin{equation}
\int_0^{\infty} e^{-(a^2+b^2) x^2}\,\cos(s x)\,\mathrm{d}x = \frac{\sqrt{\pi} e^{-s^2/4(a^2+b^2)}}{2\sqrt{a^2+b^2}} .
\end{equation}On the other hand, the integral in the second term is just the derivative of $f(s)$. Then $f(s)$ satisfies the first order inhomogenuous differential equation
\begin{equation}\label{sinedq}
f'(s)+\frac{s}{2 a^2} f(s)= \frac{b}{2 a^2 \sqrt{a^2 + b^2}} e^{-s^2/4(a^2+b^2)}.
\end{equation} Clearly the integral can be evaluated by solving this equation for $f(s)$ subject to the initial condition
\begin{equation}\label{bcondition}
f(0)=\int_0^{\infty} e^{-a^2 x^2} \mathrm{erf}(b x) \sin(0 \cdot x)\, dx = 0 .
\end{equation}
We are assured that the solution to equation \eqref{sinedq} satisfying the condition \eqref{bcondition} is the value of the given integral because $f(s)$ and its first derivative are continuous and it is known that the solution to a first order differential equation is unique provided the function and its derivative are continuous. The differential equation can be readily solved by the method of integrating factor. The method yields the general solution \begin{equation}\label{gensol}
f_g(s)=C e^{-s^2/4 a^2} + \frac{b}{2a^2\sqrt{a^2+b^2}}e^{-s^2/4a^2} \int_0^s e^{b^2 s’^2/4a^2(a^2+b^2)} ds’,
\end{equation}where $C$ is a constant of integration to be determined. Imposing the boundary condition-\eqref{bcondition} on the general solution-\eqref{gensol} yields $C=0$. By performing the change in variable $u=b s’/2a\sqrt{a^2+b^2}$, the solution assumes the integral form
\begin{equation}
f(s)= e^{-s^2/4a^2} \int_0^{b s/2\sqrt{a^2+b^2}} e^{u^2}\,\mathrm{d}u .
\end{equation}We compare the integral with the integral representation of the imaginary function,
\begin{equation}
\mathrm{erfi}(z)=\frac{2}{\sqrt{\pi}} \int_0^z e^{t^2}\,\mathrm{d}t ,
\end{equation}for all complex $z$. Then we arrive at the integral identity equation \eqref{sinintegral} for real $a$, $b$ and $c$. By analytic continuation, the result extends to complex $a$, $b$ and $c$ provided $a\neq 0$ and $\mathrm{Re}(a^2+b^2)>0$, with the square root given by its principal value.
Here is the Maple sheet for the numerical confirmation of the integral.