An algorithm is presented for calculating eigenvalues lying in the spectral gaps of the essential spectrum of Schrödinger operators. The method is applicable to the scalar and matrix Schrödinger equations. Shooting methods and Floquet theory (including a new monotonicity result for the associated Atkinson Θ-matrices) are employed to deal with the λ-dependent BVPs. The proposed method shows results with higher accuracy and lower cost comparing with those obtained using finite difference schemes combined with a contour integral λ-nonlinear eigenvalue algorithm or using the dissipative barrier scheme with domain truncations which lead to a λ-linear eigenvalue problem.