FWIR, jacobi (and SOR, Gauss-Seidel) are basically just averaging values over time steps. How many dimensions are you dealing with and what is the PDE you're solving? PDL almost surely can help you out; but even if not it'd be pretty easy to implement this in pure Perl^tm.