Some important physical effects appear due to molecular short range correlations in nematic liquid crystals. When the orientational pair potential including the higher order terms, recent theories are not yet able to give reliable calculations. This paper is devoted to developing a new numerical method to solve the equation for equili-brium state, based upon the cell theory. We do not use any empirical expansions for the orientational distribution function but obtain the accurate molecular orientational distrbutions. Then we give the accurate numerical results of the order parameter at the phase transition, the entropy change etc., and obtain the correction for the Sene-ralized mean field theory made by the short range correlations.