% GJExample1
%
digits(32); dig=digits;
ab1=vpaconvert('coeffjaclog21');
dpsi1=feval(symengine,'psi',1/4,1); dpsi1=vpa(dpsi1);
dpsi2=feval(symengine,'psi',1/2,1); dpsi2=vpa(dpsi2);
dpsi3=feval(symengine,'psi',3/4,1); dpsi3=vpa(dpsi3);
dpsi4=feval(symengine,'psi',1,1); dpsi4=vpa(dpsi4);
dpsi5=feval(symengine,'psi',5/4,1); dpsi5=vpa(dpsi5);
derr=1;
N=18;
while derr>1.5e-32
  N=N+1
  ab=ab1(1:N,:);
  xw=sgauss(dig,N,ab);
  K1=4*sum(xw(:,2)./(1+xw(:,1)));
  err=abs((K1-dpsi1+dpsi3)/(dpsi1-dpsi3));
  derr=subs(err);
end
K1=vpa(K1), derr
ab2=vpaconvert('coeffjaclog22');
derr=1;
N=18;
while derr>1.5e-32
  N=N+1
  ab=ab2(1:N,:);
  xw=sgauss(dig,N,ab);
  K2=4*sum(xw(:,2)./(1+xw(:,1)));
  err=abs((K2-dpsi2+dpsi4)/(dpsi2-dpsi4));
  derr=subs(err);
end
K2=vpa(K2), derr
ab3=vpaconvert('coeffjaclog23');
derr=1;
N=18;
while derr>1.5e-32
  N=N+1
  ab=ab3(1:N,:);
  xw=sgauss(dig,N,ab);
  K3=4*sum(xw(:,2)./(1+xw(:,1)));
  err=abs((K3-dpsi3+dpsi5)/(dpsi3-dpsi5));
  derr=subs(err);
end
K3=vpa(K3), derr