Vtyjkyuk

I am trying to approximate the digamma function in order to graph it in latex. I've already found an approximation of the gamma function from the Tex SE:
$$
Gamma(z)=
(2.506628274631sqrt(1/z) + 0.20888568(1/z)^(1.5) + 0.00870357(1/z)^(2.5) - frac(174.2106599*(1/z)^(3.5))25920 - frac(715.6423511(1/z)^(4.5))1244160)exp((-ln(1/z)-1)z)
$$
and for the digamma, $psi$, ive been using an approximation of the derivative:
$$
psi(x)=fracln Gamma(x+0.0001)-ln Gamma(x-0.0001)0.0002
$$
But im wondering if there is a better way

Yes, there is a formula
$$ψ(x) approx log(x) + frac12x - frac112x^2 +frac1120x^4 -frac1252x^6 +frac1240x^8 -frac5660x^10+frac69132760x^12 -frac112x^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) = frac1x ψ(x)$$

You can find my source here, under "Computation and Approximation".

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 = frac12$:
$$ψ(frac12) - S(frac12) = 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ψ(x+1)$$
Applying this just once for $x=frac12$, 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$.

Since you need an implementation for graphing, a simple approximation with relative low accuracy and optimized for absolute error seems suitable. A short search with computer program shows that $$psi(x) approx ln (x + a) - frac1b x, space a = 0.4849142940227510, space b=1.0271785180163817$$
approximates with an absolute error of less than $0.00123$ on $[frac12, 2^1024)$ when evaluated in IEEE-754 double precision. For arguments $x lt frac12$, the function can be computed using the reflection formula $psi(x) = psi(1-x) + pi space cot space(pi (1-x))$.