In this paper a multiharmonic balancing technique is used to develop certain algorithms to determine periodic orbits of non-linear dynamical systems with external, parametric and self excitations. Essentially, in this method the non-linear differential equations are transformed into a set of non-linear algebraic equations in terms of the Fourier coefficients of the periodic solutions which are solved by using the Newton-Raphson technique. The method is developed such that both fast Fourier transform and discrete Fourier transform algorithms can be used. It is capable of treating all types of non-linearities and higher dimensional systems. The stability of periodic orbits is investigated by obtaining the monodromy matrix. A path following algorithm based on the predictor-corrector method is also presented to enable the bifurcation analysis. The prediction is done with a cubic extrapolation technique with an arc length incrementation while the correction is done with the use of the least square minimisation technique. The under determined system of equations is solved by singular value decomposition. The suitability of the method is demonstrated by obtaining the bifurcational behaviour of rolling contact vibrations modelled by Hertz contact law.