function z=T(x,n)
y=x+sqrt(x.^2-1);
y=y.^n;
z=real(0.5*(y+1.0./y));
endfunction
