function cs = oref(s,pn,N,cs)

% cs = oref(s,pn,N,cs)
% refine location of local extrema for orthog.m

% subprograms : polydiff

% f := (1-x)^N * ( pn(x) + x^N*(1/2-x)*s((1/2-x)^2)   );
% d1 := diff(f,x);
% lprint(d1);
% d2 := diff(d1,x);
% lprint(d2);

x = cs;
pn1 = polydiff(pn);
pn2 = polydiff(pn1);
s1  = polydiff(s);
s2  = polydiff(s1);

for k = 1:8

pnx  = polyval(pn,x);
pn1x = polyval(pn1,x);
pn2x = polyval(pn2,x);
sx  = polyval(s,(1/2-x).^2);
s1x = polyval(s1,(1/2-x).^2);
s2x = polyval(s2,(1/2-x).^2);

d1 = -(1-x).^N*N./(1-x).*(pnx+x.^N.*(1/2-x).* sx ) + ...
   (1-x).^N.*(pn1x+x.^N.*N./x.*(1/2-x).*sx-x.^N.*sx+x.^N.*(1/2-x).*s1x.*(-1+2*x));

d2 = (1-x).^N.*N.^2./(1-x).^2.*(pnx+x.^N.*(1/2-x).*sx)-(1-x).^N.*N./(1-x).^2.*(pnx+x.^N.*(1/2-x).*sx)- ...
2*(1-x).^N.*N./(1-x).*(pn1x+x.^N.*N./x.*(1/2-x).*sx-x.^N.*sx+x.^N.*(1/2-x).*s1x.*(-1+2*x))+ ...
(1-x).^N.*(  ...
pn2x + x.^N.*N.^2./x.^2.*(1/2-x).*sx - x.^N.*N./x.^2.*(1/2-x).*sx - 2*x.^N.*N./x.*sx + ...
2*x.^N.*N./x.*(1/2-x).*s1x.*(-1+2*x) - 2*x.^N.*s1x.*(-1+2*x) +...
x.^N.*(1/2-x).*s2x.*(-1+2*x).^2 + 2*x.^N.*(1/2-x).*s1x ...
);

x = x - d1./d2;

end

cs = x;