How to find gain and phase margin from transfer function without Nyquist or Bode plot
Exact gain margin and phase margin finding requires us to find the roots of polynomials of \$s= 0 + j \omega\$1 constructed from the numerator and denominator polynomials of the transfer function.
Finding exact roots of high order polynomial is not possible AFAIK. But the roots can be found to very very accurate values using numerical methods. Numerical analysis software often have a roots function which can be used to find the roots numerically.
For example, to find phase margin, the equation to find the gain cross over frequency is,
$$
|N(s)| - |D(s)| = 0
$$
Where \$N, D\$ are the numerator and denominator polynomials of the open loop transfer function.
If allowed to use matlab, you can directly use the allmargin command if the transfer function or some other representation of the system is known.

Is there any easy way.

You can probably reduce the order of the system and make the polynomials lower order. But, if you are using a fast computer, it may not be worth the effort.
1 \$s = \sigma + j \omega\$ and lies on the typical s plane contour.