function y=f2(x,m,beta,c)
if m==0
  p=1;
elseif m==1
  p=x-c;
else
  if m==2*floor(m/2), error('invalid m'), end
  p=ones(size(x,1),1);
  for k=1:(m-1)/2
    p=(x.^2+((2*k-1)*pi/beta)^2).*p;
  end
  p=(x-c).*p;
end
y=p.*sinh(beta*x)./((1+cosh(beta*x)).*(x-c));
