The cylindrical Bessel differential equation and the spherical Bessel differential equation in the interval with Neumann boundary conditions are considered. The eigenfunctions are linear combinations of the Bessel function or linear combinations of the spherical Bessel functions. The orthogonality relations with analytical expressions for the normalization constant are given. Explicit expressions for the Lommel integrals in terms of Lommel functions are derived. The cross product zeros and are considered in the complex plane for real as well as complex values of the index and approximations for the exceptional zero are obtained. A numerical scheme based on the discretization of the two-dimensional and three-dimensional Laplace operator with Neumann boundary conditions is presented. Explicit representations of the radial part of the Laplace operator in form of a tridiagonal matrix allow the simple computation of the cross product zeros.