A quasi-Newton method with cubic regularization is designed for solving Riemannian unconstrained nonconvex optimization problems. The proposed algorithm is fully adaptive with at most ${\cal O} (\epsilon_g^{-3/2})$ iterations to achieve a gradient smaller than $\epsilon_g$ for given $\epsilon_g$, and at most $\mathcal O(\max\{ \epsilon_g^{-\frac{3}{2}}, \epsilon_H^{-3} \})$ iterations to reach a second-order stationary point respectively. Notably, the proposed algorithm remains applicable even in cases of the gradient and Hessian of the objective function unknown. Numerical experiments are performed with gradient and Hessian being approximated by forward finite-differences to illustrate the theoretical results and numerical comparison.