The Bessel function(s) are pretty difficult to calculate and generally have to be calculated numerically via any of several methods. Looking at the numbers you posted, I'd say that Jmat is probably using a different algorithm and/or precision level than the others. My hunch would be that it's using a series expansion, while the other software is using one of the more clever algorithms.

PS: In Mathematica, you should probably use "N[..., 20]" instead of "SetPrecision". It not quite the same and using SetPrecision can cause odd effects if you're not careful, like accidentally triggering numerical evaluation too early.

It looks like Jmat is using the series expansion (about zero) for Re(z)^2 < (Re(n)+1)/10, with z the argument and n the order. They switch to the asymptotic form for larger arguments, but they're only using the first term.

To improve accuracy, they could truncate the large order expansion when the terms start growing or drop below precision. See the wonderful DLMF S10.17 (http://dlmf.nist.gov/10.17) for the higher order terms.

Note that this expansion will not perform well for large order, and they should switch to one of the series in order instead. Olver's uniform expansion is lovely, but a bit tedious to implement.

Ok, improving the Bessel functions is now the first priority thing to fix :)