The error function is defined as
$$\operatorname{erf}(x)=\frac 2 {\sqrt \pi}\int_{0}^x e^{-t^2}dt$$
It is not an elementary function. Since from the definition it is immediate (FTCI) that
$$\operatorname{erf}'(x)=\frac 2 {\sqrt \pi}e^{-x^2}$$
the primitive of $e^{-x^2}$ is expressible as
$$\int e^{-x^2} dx =\frac{\sqrt \pi}{2}\operatorname{erf}(x)-C $$
since any two primitives of a function $f$ differ by a constant (FTCII)
As a consequence your primitive can't be expressed in terms of elementary functions.
$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
With $\ds{r \equiv {1 + \root{3}\ic \over 2} = \expo{\pi\ic/3}}$
\begin{align}&\color{#c00000}{
\int_{0}^{1}{\ln\pars{\ln\pars{1/x}} \over x^{2} - x + 1}\,\dd x}
=\int_{0}^{1}{\ln\pars{\ln\pars{1/x}} \over \pars{x - r}\pars{x - r^{*}}}\,\dd x
\\[3mm]&=\int_{0}^{1}\ln\pars{\ln\pars{1/x}}
\pars{{1 \over x - r} - {1 \over x - r^{*}}}\,{1 \over r - r^{*}}\,\dd x
\\[3mm] & =
{1 \over \Im\pars{r}}\,\Im\int_{0}^{1}
{\ln\pars{\ln\pars{1/x}} \over x - r}\,\dd x
\end{align}
With $\ds{x \equiv \expo{-t}}$:
\begin{align}&\color{#c00000}{
\int_{0}^{1}{\ln\pars{\ln\pars{1/x}} \over x^{2} - x + 1}\,\dd x}
={2\root{3} \over 3}\Im\int_{\infty}^{0}{\ln\pars{t} \over \expo{-t} - r}
\,\pars{-\expo{-t}\,\dd t}
\\[3mm]&=-\,{2\root{3} \over 3}\Im\bracks{{1 \over r}\int_{0}^{\infty}
{\ln\pars{t}\expo{-t} \over 1 - \expo{-t}/r}\,\dd t}
\\[3mm]&=-\,{2\root{3} \over 3}\Im\bracks{{1 \over r}
\sum_{n = 1}^{\infty}{1 \over r^{n - 1}}\int_{0}^{\infty}
\ln\pars{t}\expo{-nt}\,\dd t}
\\[3mm]&=-\,{2\root{3} \over 3}\Im\lim_{\mu\ \to\ 0}\partiald{}{\mu}\bracks{
\sum_{n = 1}^{\infty}r^{-n}\int_{0}^{\infty}t^{\mu}\expo{-nt}\,\dd t}
\\[3mm]&=-\,{2\root{3} \over 3}\Im\lim_{\mu\ \to\ 0}\partiald{}{\mu}\bracks{
\sum_{n = 1}^{\infty}{r^{-n} \over n^{\mu + 1}}\Gamma\pars{\mu + 1}}
\\[3mm]&=-\,{2\root{3} \over 3}\Im\lim_{\mu\ \to\ 0}\partiald{}{\mu}\bracks{
\Gamma\pars{\mu + 1}{\rm Li}_{\mu + 1}\pars{r^{*}}}
\\[3mm]&=-\,{2\root{3} \over 3}\Im\lim_{\mu\ \to\ 0}\partiald{}{\mu}\bracks{
\Gamma\pars{\mu + 1}{\rm Li}_{\mu + 1}\pars{\expo{-\pi\ic/3}}}
\end{align}
$\ds{{\rm Li}_{1}\pars{z} = -\ln\pars{1 - z}}$. Derivatives of the
PolyLogarithm, respect of the order, can be evaluated from its integral representation.
Also, see Hurwitz Zeta Function.
Best Answer
Presented below is a full solution via elementary integration.
The polynomial $x^4-2x^3-6x-1$ admits two real roots, $a=$ -0.1650 and $b=$ 2.8068, computed with the standard quartic-root formulae below
$$a,b=\frac12\left(1+\sqrt r \pm \sqrt{3-r+{14}{r^{-\frac12}}}\right)\tag1 $$ with $ r= 1+ \left(\frac43\right)^{\frac23}\left( \sqrt[3]{\sqrt{129}+9} - \sqrt[3]{\sqrt{129}-9} \right) $. Thus, the denominator of the integrand factorizes as
$$x^4-2x^3-6x-1=(x-a)(x-b)[x^2+(a+b-2)x-(ab)^{-1}] $$ leading to following partial fractionalization \begin{align} \frac{(x^2-1)(x^2+3)}{x^4-2x^3-6x-1} =1+ \frac {p}{x-a} + \frac {q}{x-b}+ \frac{(2-p-q )x-\frac1{ab}(2+\frac {p}a+ \frac {q}b)}{x^2+(a+b-2)x-\frac1{ab}} \end{align} with $p=\frac{a^3+a^2+3a-1}{2a^3-3a^2-3}$ and $q=\frac{b^3+b^2+3b-1}{2b^3-3b^2-3}$. Then, integrate to obtain
\begin{align} & \int \frac{(x^2-1)(x^2+3)}{x^4-2x^3-6x-1}dx\\ =& \>x +2p \ln|x-a| +2q \ln|x-b|+(1-p-q)\ln\left(x^2+(a+b-2)x-\frac1{ab}\right)\\ & \>-\frac{(1-p-q)(a+b-2)+\frac2{ab}(1+\frac pa + \frac qb)}{\sqrt{ -\frac{(a+b-2)^2}4-\frac1{ab}}} \tan^{-1}\frac{x+\frac{a+b-2}2}{\sqrt{ -\frac{(a+b-2)^2}4-\frac1{ab}}}+C \end{align} The result, though rather involved in appearance, is expressed in terms of $a$ and $b$, the two real roots given in (1).