function y=f1(x,m,beta)
if m==0
  p=1;
else
  p=ones(size(x,1),1);
  for k=1:m/2
    p=(x.^2+((2*k-1)*pi/beta)^2).*p;
  end
end
y=p.*sinh(beta*x)./(1+cosh(beta*x));
