There is no hope as such for Mathematica problem. I would have to reply solely on IDL and probable integration of which with Mathematica. I still have one more option open, to write iteration method inside Mathematica. One such fine example can be found here:
pi[n_, theta_] := 
     p1=p2*(2i-1)/(i-1) Cos[theta] –  i/(i-1) p1;

But it doesn’t make a sense to try too much with Mathematica, when IDL is working pretty good.