As a reminder for everybody, the general Levin transformation for a sequence of partial sums $S_n=\sum\limits_{j=0}^n a_j$ looks like this:
$$\mathcal{L}_k^{(n)}=\frac{\Delta^k\left((n+b)^{k-1}\frac{S_n}{g(n)}\right)}{\Delta^k\left((n+b)^{k-1}\frac1{g(n)}\right)}$$
where $\Delta^k$ is the usual $k$-th forward difference operator, $g(n)$ is an auxiliary sequence, and $b$ is an adjustable real parameter that is not a nonpositive integer; more explicitly, we have
$$\mathcal{L}_k^{(n)}=\frac{\sum\limits_{j=0}^k (-1)^j\binom{k}{j}(n+j+b)^{k-1}\frac{S_{n+j}}{g(n+j)}}{\sum\limits_{j=0}^k (-1)^j\binom{k}{j}(n+j+b)^{k-1}\frac1{g(n+j)}}$$
Often, $b$ is taken to be $1$, and $n$ is taken to be $0$. The various Levin transformations correspond to different choices of the auxiliary sequence $g(n)$; to give two of the simplest cases, the $u$-transformation, for instance, takes $g(n)=(n+b)a_n$, while the $t$-transformation uses $g(n)=a_n$.
The article you should be looking at, apart from David Levin's original paper, is E.J. Weniger's Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series; in there, he gives a short FORTRAN routine for implementing Levin's transformations. A more elaborate implementation by Fessler, Ford, and Smith is available at Netlib.
Just to show how easy it is to implement the Levin transformation, here's a demonstration program I wrote for the TI-83 Plus calculator for summing the alternating harmonic series $\sum\limits_{k=1}^\infty \frac{(-1)^{k+1}}{k}$, based on Weniger's short FORTRAN routine (comments are delimited by a backslash):
PROGRAM:XTRPOL
Prompt N \\ number of terms of the series to use
1→P:0→S \\ P: alternating sign S: stores partial sums
For(K,1,N)
P/K→U \\ K-th term of the series, change to sum a different series
S+U→S
1→B \\ adjustable parameter for the Levin transformation
(B+K-1)U→T \\ Levin u-transform; for t version, remove the (B+K-1)
prgmLEVINT
Pause Y
-P→P
End
Y
which uses the subroutine
PROGRAM:LEVINT
(B+K-1)ֿ¹→W
W/T→U
If K=1
ClrList ∟DL,∟NL
U→∟DL(K)
SU→∟NL(K)
1-W→V
For(J,K-1,1,-1)
(B+J-1)W→U
∟NL(J+1)-U∟NL(J)→∟NL(J)
∟DL(J+1)-U∟DL(J)→∟DL(J)
WV→W
End
10^(99)→Y
If abs(∟DL(1))≥10^-99
∟NL(1)/∟DL(1)→Y
Y
I haven't bothered to implement the stopping rules described by Fessler, Ford, and Smith in that demo program, but it's doable. Translating that short routine to your language of choice should be straightforward.
As a note, the algorithm here looks very simple, due to the exploitation of the recursive identities satisfied by the forward differences.
