I'm trying to efficiently graph the function $I(t)$ for $t>0$ where
$$I(t) :=\int_0^t e^{-\lambda s} \,\text{erf}(\ln(t-s))\,\text{d}s,\qquad \lambda>0$$
but its evaluation is beyond my powers of integration. I can use a numerical integration package to plot it, but it is pretty slow. Happy to compute a decent approximation if it's available.
Mathematica was not able to provide an answer, and I've tried some integral substitutions like $u=\ln(t-s)$, as well as looked through the fantastic table of erf integrals, but no solution seems to be clear.
My current strategy is to approximate the exponential by a quartic polynomial $$\exp(-\lambda x)\approx \sum_{k=0}^4 c_k (\lambda x)^k$$ in some region $x\in[0,R]$, and set to zero when $x>R$. Then I can compute $\int s^k \text{erf}(\ln(t-s))\,\text{d}s$, but there are a painful number of terms when using a quartic - are there any better suggestions? Perhaps some other approximation of the $\exp$, or nice series expansion of the function?
Thanks!
Edit: Just to be clear, imagine I am given ten thousand different values of $\lambda$, and I want to graph $I(t)$ for each of them in the region $t\in [0,100]$. How can I do this efficiently?