Yes, there is a formula
$$ψ(x) \approx \log(x) - \frac{1}{2x} - \frac{1}{12x^2} +\frac{1}{120x^4} -\frac{1}{252x^6} +\frac{1}{240x^8} -\frac{5}{660x^{10}}+\frac{691}{32760x^{12}} -\frac{1}{12x^{14}}$$
This is especially accurate for larger values of $x$. If $x$ is small, you can shift $x$ to a higher value using the relation
$$ψ(x+1) = \frac{1}{x} + ψ(x)$$
You can find my source here, under the "Computation and Approximation" section.
Edit: Accuracy of Approximation
This series is quite accurate for $x$ in the interval $(1,\infty)$. Given the OP's bounds of $x = (0,20)$, we can find the difference between the expected and calculated results (note that $S(x)$ is the series in question)
$$ψ(1) - S(1) = 0.0674022$$
$$ψ(20) - S(20) = 4.44089×10^{-16}$$
From these results we find that we get exactly what we expected: the series is decent in this interval, but is MUCH better for larger $x$. However, we get a problem for $x = \frac{1}{2}$:
$$ψ(\frac{1}{2}) - S(\frac{1}{2}) = 1285.81$$
Obviously we are way off here. However, some context is useful. Our series blows up for $x<1$ towards $-\infty$, and so we can't expect anything useful here. Nevertheless, we can use the shift formula above to mitigate this issue, rewriting it as
$$ψ(x) = ψ(x+1) - \frac{1}{x}$$
Applying this just once for $x=\frac{1}{2}$, we get a difference of $1.9816604...$, an approximation $\approx 650$ times better. Thus, we see that we can use the series provided to get a decent approximation of the graph of $ψ(x)$ for all $x>0$.
Quoting from Wikipedia: In mathematics, an elementary function is a function of a single variable (typically real or complex) that is defined as taking sums, products, and compositions of finitely many polynomial, rational, trigonometric, hyperbolic, and exponential functions, including possibly their inverse functions (e.g., $\arcsin$, $\log$, or $x^{1/n}$).
As you can see, the exponential function is on the list. Now, the modified Bessel function $I_{0}$ is a non-elementary function. You can express it in terms of elementary functions but such an expression will always involve an infinite process, e.g., infinite series, integration, etc. It may be hard to accept that adding an extra factorial to the sum changes things so dramatically, but that is the case here. You have to accept that what you are asking for is not possible, not because mathematicians are incompetent but because of the nature of mathematics.
It was shown by Hankel that
$$
I_0 (2\sqrt x ) \sim\frac{{\mathrm{e}^{2\sqrt x } }}{{2\sqrt{\pi} x^{1/4} }}\left( {1 + \frac{1}{{16\sqrt x }} + \frac{9}{{512x}} + \cdots } \right),
$$
as $x\to +\infty$. This is a divergent asymptotic expansion in the sense of Poincaré. For the full expansion, see $(10.40.1)$ in the NIST DLMF. The series is divergent for all $x>0$ but any truncated version of it yields a good approximation for $I_0$ for large positive values of $x$. This asymptotics also shows you how the function deviates from $\mathrm{e}^x$.
When $x<0$, your sum can be written in terms of the $J_0$ Bessel function as $J_0 (2\sqrt { - x} )$. The behaviour for large negative $x$ is given by another result of Hankel and it shows an oscillatory behaviour. At leading order
$$
J_0 (2\sqrt { - x} ) \approx \frac{1}{{\sqrt \pi ( - x)^{1/4} }}\cos \left( {2\sqrt { - x} - \frac{\pi }{4}} \right),
$$
when $x$ is large and negative.
The above two asymptotics shows that the behaviour of your series is rather different when $x>0$ or $x<0$ and may also convince you that it is not possible to capture two such different behaviours by a single approximation.
Best Answer
The Digamma function is analytic at $1$, with Taylor series
$$ \psi(z) = -\gamma - \sum_{n=1}^\infty \zeta(n+1) (1-z)^n $$
This converges for $|z-1|<1$.
EDIT: $\psi$ is analytic in $\mathbb C$ except for poles at the nonpositive integers. So $\psi(1/z)$ is analytic in the disk of radius $1$ around $1$, and its Taylor series around $1$ converges there. So if this Taylor series is $$\psi(1/z) = \sum_{n=0}^\infty c_n (z-1)^n$$, we can write $$ \psi(x) = \sum_{n=0}^\infty c_n \left(\frac{1}{x}-1\right)^n $$ and use partial sums of this as our approximations. For example, the sum up to $n=6$ is
$$-\gamma -\frac{\pi^{2} \left(\frac{1}{x}-1\right)}{6}+\left(\frac{\pi^{2}}{6}-\zeta \! \left(3\right)\right) \left(\frac{1}{x}-1\right)^{2}+\left(-\frac{\pi^{2}}{6}+2 \zeta \! \left(3\right)-\frac{\pi^{4}}{90}\right) \left(\frac{1}{x}-1\right)^{3}+\left(\frac{\pi^{2}}{6}-3 \zeta \! \left(3\right)+\frac{\pi^{4}}{30}-\zeta \! \left(5\right)\right) \left(\frac{1}{x}-1\right)^{4}+\left(-\frac{\pi^{2}}{6}+4 \zeta \! \left(3\right)-\frac{\pi^{4}}{15}+4 \zeta \! \left(5\right)-\frac{\pi^{6}}{945}\right) \left(\frac{1}{x}-1\right)^{5}+\left(\frac{\pi^{2}}{6}-5 \zeta \! \left(3\right)+\frac{\pi^{4}}{9}-10 \zeta \! \left(5\right)+\frac{\pi^{6}}{189}-\zeta \! \left(7\right)\right) \left(\frac{1}{x}-1\right)^{6} $$ where the error is about $0.0020$ at $x=2$ and $0.0206$ at $x=3$.