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_] := 
Module[{p1,p2},
   p1=0;
   p2=1;
   Do[
     p1=p2*(2i-1)/(i-1) Cos[theta] –  i/(i-1) p1;
     tmp=p2;
     p2=p1;
     p1=tmp,{i,2,n}]
   p2
]

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