The relevant lines of code are like this. The last calculation takes 1 hour 20 min. >LegendreQmy[np_,xp_]:=LegendreP[np,xp]*Log[(xp+1)/(xp1)]/2Sum[LegendreP[ksum,xp]*LegendreP[npksum1,xp]/(npksum),{ksum,0,np1}]/;np>1 >IntPP[kp_,lp_,xp_]:=2*Sum[(App[(kp+lp+mp)/2kp]*App[(kp+lp+mp)/2lp]*App[(kp+lp+mp)/2mp]/App[(kp+lp+mp)/2])*((2*mp+1)/(kp+lp+mp+1))*LegendreQmy[mp,xp],{mp,Abs[kplp],kp+lp,2}] >Clear ;b[gp_,np_]:=(2*np+1)*(gp^np)
>Clear[zet,etap,phip];zet[etap_,phip_]:=1+(1etap)*Exp
>Clear[gp,zp,RplusmatrixLeg];RplusmatrixLeg[zp_,gp_]=Table[Sqrt[b[gp,ik1]*b[gp,jk1]]*IntPP[ik1,jk1,zp],{ik,Nc},{jk,Nc}];
>g=0.8;etad=0.01;dphi=2*Pi/90;Clear;Rleg=Table[Rini,{ietad,1,101},{iphi,1,46}];Do[Do[Rleg[[ietad,iphi]]=RplusmatrixLeg[zet[etad*(ietad1),Pi+dphi*(iphi1)],g],{iphi,1,46}],{ietad,1,100}]; Do[Rleg[[ietad,1]]=RplusmatrixLeg[zet[etad*(ietad1),Pi+0.000001],g],{ietad,1,100}]
