function F=pFe(z,q)
N=100; H=zeros(N,N); V=zeros(N,1);
for m=1:N-1
  H(m,m)=(2*m-1)^2; H(m,m+1)=q; H(m+1,m)=q;
end
H(N,N)=(2*N-1)^2; H(1,1)=H(1,1)+q; E=eig(H); a=E(1);
V(1)=1; V(2)=(a-1-q)/q;
for m=2:N-1
  V(m+1)=(a-(2*m-1)^2)/q*V(m)-V(m-1);
end
[C,m]=min(abs(V)); A=V(1:m)/norm(V(1:m),'fro');
F=0;
for r=1:m
  F=F+(2*r-1)*A(r)*bessely(2*r-1,2*sqrt(q)*sinh(z));
end
F=ce1(eps,q)/(sqrt(q)*A(1))*coth(z).*F;