In this paper, we present a detailed numerical study of the non-linear dynamics in structural components under unilateral contact constraints. Here, the unilateral term characterises the constitutive law of the restoring force in the constraints as they only sustain elastic reactions in one direction, either compressive or tensile. Thus, the non-differentiability of the contact law at the discontinuity point is the only source of non-linearity. In our approach, the discrete lattice method (DLM) is used to treat the continuous system as a piecewise linear model. Thus, the trajectory of each node in the discrete model would be a sequence of smooth solutions with the switching times between them. The application of the one-step integration scheme allows us to detect the occurrence of contact (i.e. the instants that the lattice nodes cross the discontinuity boundary) and consequently update the active constraints. We also consider embedding the bisection algorithm into the time integration procedure to localise the instants at which the nodes cross the boundary and minimise the accumulative error. Subsequently, the resulting unconditionally stable integration scheme is utilised as the modelling tool in combination with the shooting technique to perform a novel non-smooth modal analysis. In analogy with the smooth non-linear systems, the evolution of non-smooth periodic motions is presented in the frequency-stiffness plots. We apply our method to obtain non-linear normal modes (NNMs) for a number of representative problems, including a bar-obstacle system, a beam-substrate system and a granular chain with tensionless interactions. These numerical examples demonstrate the efficiency of the solution procedure to trace the family of energy-independent non-linear modes across the range of contact stiffnesses. Moreover, the stability analysis of the modes on the plot backbone reveal that they may become unstable due to the interaction with the higher modes or bifurcation of a new NNM with a period equal to the integer multiple of the main mode period. Our results also indicate that the bilinear formula can accurately predict the non-linear frequencies only if the corresponding mode exhibits a smooth character, regardless of the commutativity conditions of the system stiffness matrix. However, it is obvious that the assumption of smooth bilinear behaviour for non-linear modes is not generally valid. This highlights the importance of the present numerical framework for the computation of non-smooth resonance frequencies.