Ostrogradsky-Hermite method. As I have written in more detail here, we can use the Ostrogradski-Hermite method to integrate a rational fraction $P(x)/Q(x)$ without decomposing it into partial fractions and without finding the multiple roots of the denominator. You can also find a simple description of this method here, subsection 4.5.2, as well as some exercises (1301-1304).
Assume that $\deg P<\deg Q$. Then, there exist polynomials $$P_{1},\quad Q_{1}=\gcd \left\{ Q,Q^{\prime }\right\} ,\quad P_{2}\quad \text{and}\quad Q_{2}=Q/Q_{1},$$ with $\deg P_{1}<\deg Q_{1}$, $\deg P_{2}<\deg Q_{2}$, such that
\begin{equation*}
\int \frac{P}{Q}\, dx=\frac{P_{1}}{Q_{1}}+\int \frac{P_{2}}{
Q_{2}}\, dx.\tag{1}
\end{equation*}
Since the integrand is
\begin{equation*}
\frac{P}{Q}=\frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left(
x^{3}-x+1\right) ^{2}},\tag{2}
\end{equation*}
and
\begin{eqnarray*}
Q &=&\left( x^{3}-x+1\right) ^{2} \\
Q^{\prime } &=&2\left( x^{3}-x+1\right) \left( 3x^{2}-1\right) \\
Q_{1} &=&\gcd \left\{ Q,Q^{\prime }\right\} =x^{3}-x+1,\\
Q_{2}&=&\frac{Q}{Q_{1}}=x^{3}-x+1,
\end{eqnarray*}
we can expand the given integral as
\begin{equation*}
\int \frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left( x^{3}-x+1\right) ^{2}}\, dx=
\frac{Ax^{2}+Bx+C}{x^{3}-x+1}+\int \frac{Dx^{2}+Ex+F}{x^{3}-x+1}\, dx,\tag{3}
\end{equation*}
where $A$, $B$, $C$, $D$, $E$, $F$ are constants. They can be found by
differentiating $(3)$, reducing both sides to a common denominator and equating the coefficients of like powers of $x$ of the numerators. We thus have
\begin{eqnarray*}
\frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left( x^{3}-x+1\right) ^{2}} &=&
\frac{\left( 2Ax+B\right) \left( x^{3}-x+1\right) -\left( 3x^{2}-1\right)
\left( Ax^{2}+Bx+C\right) }{\left( x^{3}-x+1\right) ^{2}}\\ &&+\frac{\left(Dx^{2}+Ex+F\right)\left(x^{3}-x+1\right)}{\left( x^{3}-x+1\right) ^{2}}.
\end{eqnarray*}
Since $\left( x^{2}-3x+1/3\right) ^{2}=x^{4}-6x^{3}+(29/3)x^{2}-2x+1/9$, after some algebra we obtain
\begin{align}
x^{4}-6x^{3}+(29/3)x^{2}-2x+1/9&= Dx^{5}+\left( -A+E\right) x^{4}+\left( -2B-D+F\right) x^{3}\\
&\qquad+\left(-A-3C+D-E\right) x^{2}+\left( 2A+E-F\right) x\\
&\qquad+\left( B+C+F\right),
\end{align}
which means that the constants satisfy the following system of equations
\begin{equation*}
\left\{
\begin{array}{l}
D =0,\\ -A+E=1,\\ -2B-D+F=-6, \\
-A-3C+D-E =\frac{29}{3},\\ 2A+E-F=-2, \\ B+C+F=\frac{1}{9},
\end{array}
\right.
\end{equation*}
whose solution is
\begin{equation*}
A=-1,\, B=3,\, C=-\frac{26}{9},\, D=0,\, E=0,\, F=0.
\end{equation*}
Hence
\begin{equation*}
\int \frac{\left( x^{2}-3x+1/3\right) ^{2}}{\left( x^{3}-x+1\right) ^{2}}\, dx=
\frac{-x^{2}+3x- 26/9 }{x^{3}-x+1}+C.\tag{4}
\end{equation*}