Conventional nonlinear beam elements often fall short in capturing the full range of torsional, axial, and lateral deformations, particularly in structures experiencing complex helical buckling. In this work, we extend an existing finite element formulation [1] to investigate the twisting and buckling behaviour of slender rods with nonuniform cross-sections. Our method is based on a modified Euler–Bernoulli beam theory. Symbolic derivation and automated element assembly were performed in Mathematica, followed by extensive parametric studies in MATLAB. To assess the reliability and robustness of the model, we conducted a mesh sensitivity analysis, comparing our custom MATLAB implementation with results from the commercial finite element software Abaqus. A parametric analysis further explored how a geometric parameter, such as cross-sectional radius, affects buckling behaviour. The results reveal consistent patterns and strong agreement between our model and independent simulations. Overall, the approach is both computationally efficient and reliable, making it a strong candidate for integration into broader finite element libraries, with potential applications in soft robotics, coiled tubing, and biomechanical systems.