In this paper, an elegant and easy to implement numerical method using matrix mechanics approach is proposed, to solve the time independent Schrodinger equation (TISE) for Morse potential. It is specifically applied to non-homogeneous diatomic molecule HCl to obtain its rotating-vibrator spectrum. While matrix diagonalization technique is utilised for solving TISE, model parameters for Morse potential are optimized using variational Monte-Carlo (VMC) approach by minimizing χ 2 − value. Thus, validation with experimental vibrational frequencies is completely numerical based with no recourse to analytical solutions. The ro-vibrational spectra of HCl molecule obtained using the optimized parameters through VMC have resulted in least χ 2 − value as compared to those determined using best parameters from multiple regression analysis of analytical expressions. Numerical algorithm for solving the Hamiltonian matrix has been implemented utilizing Free Open Source Software (FOSS) Scilab and simulation results are matching well with those obtained using analytical solutions from Nikiforov-Uvarov (NU) method and asymptotic iteration method (AIM).