We develop a new Finite Element to accurately model twisting of rods and capture the bifurcation scenarios leading to helical buckling and various further post-buckling states. Since standard nonlinear beam elements do not account for nonlinearities in torsional modes, as well as for coupling between axial, lateral and torsional modes, we derive a new beam element, which allows us to describe complex helical buckling bifurcation scenarios of a rod subjected to a twisting load. The formulated beam element is systematically tested to assess its predictive capabilities in determining critical torsional buckling loads and its sensitivity to a number of elements used. Once the model is validated against commercial FE software (ABAQUS), we focus our attention on computing bifurcation scenarios to observe various complex helical configurations and transitions between them. The analysis reveals coexistence between helices with multiple loops for certain values of twisting load. Additionally, we trace the transition onsets between stable helical configurations. The developed FE can be applied to study complex buckling mechanics of engineering and biological structures.