Abstract: Strongly interacting particles (hadrons) consist of baryons and mesons. The strong interaction is a very complicated problem, however, it simplifies greatly in the sector of heavy quarkonium. The first heavy quarkonium states discovered in the early seventies of the last century were the charmonium states consisting of c quark and ̅ antiquark. Later the bottomonium states consisting of b quark and ̅ antiquark were discovered. In the last decade of the last century, the mesons consisting of b quark and ̅ antiquark was also discovered. These states of heavy quarkonium helped us understand many features of the theory of quantum Chromodynamics (QCD) which describe the interactions between quarks and gluons. However, since the theory of QCD proved to be very difficult to solve, many models were proposed to explain the spectrum and other properties of the heavy quarkonium states In this thesis we have used matrix methods to solve a quarkonium model with Hulthen potential describing the short range interaction. We evaluate the spectrum of heavy quarkonium by solving the non-relativistic Schrodinger equation for various orbital angular momentum quantum numbers. We take spin-spin, spin-orbit and tensor interactions into account. We have compared our results with experimental results and also with the results of other relevant models. We conclude that this model gives a good description of the spectrum of heavy quarkonium states.