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)/(xp-1)]/2-Sum[LegendreP[ksum,xp]*LegendreP[np-ksum-1,xp]/(np-ksum),{ksum,0,np-1}]/;np>1 >IntPP[kp_,lp_,xp_]:=2*Sum[(App[(kp+lp+mp)/2-kp]*App[(kp+lp+mp)/2-lp]*App[(kp+lp+mp)/2-mp]/App[(kp+lp+mp)/2])*((2*mp+1)/(kp+lp+mp+1))*LegendreQmy[mp,xp],{mp,Abs[kp-lp],kp+lp,2}] >Clear ;b[gp_,np_]:=(2*np+1)*(gp^np)
>Clear[zet,etap,phip];zet[etap_,phip_]:=1+(1-etap)*Exp
>Clear[gp,zp,RplusmatrixLeg];RplusmatrixLeg[zp_,gp_]=Table[Sqrt[b[gp,ik-1]*b[gp,jk-1]]*IntPP[ik-1,jk-1,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*(ietad-1),-Pi+dphi*(iphi-1)],g],{iphi,1,46}],{ietad,1,100}]; Do[Rleg[[ietad,1]]=RplusmatrixLeg[zet[etad*(ietad-1),-Pi+0.000001],g],{ietad,1,100}]
|