%% Author: M. Turhan ÇOBAN
%% Created: 2021-11-19

function [x,C] = root_P(n)
 x=zeros(n,1);
 C=zeros(n,1);
 for k=1:n
   %% initial root finding formula
   theta_n_k=(n-k+3.0/4.0)/(n+0.5)*pi;
   sintk=sin(theta_n_k);
   costk=cos(theta_n_k);
   xk=(1.0-1.0/(8.0*n*n)+1.0/(8.0*n*n*n)-1.0/(384.0*n*n*n*n)*(39.0-28.0/(sintk*sintk)))*costk;
   dnu=1.0;
   eps=1.0e-10;
   while(dnu>eps)
     %% Newton-Raphson root finding
     Pn=P1(n,xk); 
     Pn_1=P1(n-1,xk);
     xold=xk;
     dsk=(1.0-xk*xk)/n*Pn/(Pn_1-xk*Pn);
     xk=xk-dsk;
     dnu=abs(xold-xk);
   end
   x(k)=xk;
   Pn=P1(n,xk);
   Pn1=Pn_1*n;
   C(k)=2.0*(1.0-xk*xk)/(Pn1*Pn1);
 endfor
endfunction