% FHL A function used in Example3.64.
%
function [y,N]=fHL(t,x,eps0)
squ=sqrt(2*x*t); y0=sinh(squ)*sin(squ);
y=y0; y1=0; N=0;
while abs((y1-y)/y)>eps0
  y1=y; N=N+10;
  th=linspace(0,pi,N+1); th1=th(2:N);
  y=(y0+sum(exp(squ*cos(th1)).*sin(squ*cos(th1))))/N;
end
