Generate the contour plots of the modulus and phase of the Hankel function $$H_0^{(1)}(z)$$ [1].

Create a grid of values for the domain.

[X,Y] = meshgrid(-4:0.002:2,-1.5:0.002:1.5);

Calculate the Hankel function over this domain and generate the modulus contour plot.

H = besselh(0,X+1i*Y);
contour(X,Y,abs(H),0:0.2:3.2)
hold on

In the same figure, add the contour plot of the phase.

contour(X,Y,rad2deg(angle(H)),-180:10:180)
hold off

Plot the real and imaginary parts of the Hankel function of the second kind and examine their asymptotic behavior.

Calculate the Hankel function of the second kind ${\mathit{H}}_0^{\left(2\right)}\left(\mathit{z}\right)={\mathit{J}}_0\left(\mathit{z}\right)-\mathit{i}{\mathit{Y}}_0\left(\mathit{z}\right)$ in the interval $\left[0.1,25\right]$.

k = 2;
nu = 0;
z = linspace(0.1,25,200);
H = besselh(nu,k,z);

Plot the real and imaginary parts of the function. In the same figure, plot the linear combination $\sqrt{{\mathit{J}}_0^2\left(\mathit{z}\right)+{\mathit{Y}}_0^2\left(\mathit{z}\right)}$, which reveals the asymptotic behavior of the magnitudes of the real and imaginary parts.

plot(z,real(H),z,imag(H))
grid on
hold on
M = sqrt(real(H).^2 + imag(H).^2);
plot(z,M,'--')
legend('$J_0(z)$', '$Y_0(z)$', '$\sqrt{J_0^2 (z) + Y_0^2 (z)}$','interpreter','latex')

Calculate the exponentially scaled Hankel function ${\mathit{H}}_1^{\left(2\right)}\left(\mathit{z}\right){\cdot \mathit{e}}^{\mathrm{iz}}$ on the complex plane and compare it to the unscaled function.

Calculate the unscaled Hankel function of the second order on the complex plane. When z has a large positive imaginary part, the value of the function quickly diverges. This phenomenon limits the range of computable values.

k = 2;
nu = 1;
x = -5:0.4:15;
y = x';
z = x + 1i*y;
scaled = 1;
H = besselh(nu,k,z);
surf(x,y,imag(H))
xlabel('real(z)')
ylabel('imag(z)')

Now, calculate ${\mathit{H}}_1^{\left(2\right)}\left(\mathit{z}\right){\cdot \mathit{e}}^{\mathrm{iz}}$ on the complex plane and compare it to the unscaled function. The scaled function increases the range of computable values by avoiding overflow and loss of accuracy when z has a large positive imaginary part.

Hs = besselh(nu,k,z,scaled);
surf(x,y,imag(Hs))
xlabel('real(z)')
ylabel('imag(z)')