A method to compute the numerical derivative of eigenvalues of parameterized crystal field Hamiltonian matrix is given, based on the numerical derivatives the general iteration methods such as Levenberg-Marquardt, Newton method, and so on, can be used to solve crystal field parameters by fitting to experimental energy levels. With the numerical eigenvalue derivative, a detailed iteration algorithm to compute crystal field parameters by fitting experimental energy levels has also been described. This method is used to compute the crystal parameters of Yb^3＋ in Sc2O3 crystal, which is prepared by a co-precipitation method and whose structure was refined by Rietveld method. By fitting on the parameters of a simple overlap model of crystal field, the results show that the new method can fit the crystal field energy splitting with fast convergence and good stability.
