From the equation of motion of a "bare" non-uniform beam (without any concentrated elements), an eigenfunction in term of four unknown integration constants can be obtained. When the last eigenfunction is substituted into the three compatible equations, one force-equilibrium equation, one governing equation for each attaching point of the concentrated element, and the boundary equations for the two ends of the beam, a matrix equation of the form [B]{C} = {0} is obtained. The solution of |B| = 0 (where ${mid}{cdot}{mid}$ denotes a determinant) will give the "exact" natural frequencies of the "constrained" beam (carrying any number of point masses or/and concentrated springs) and the substitution of each corresponding values of {C} into the associated eigenfunction for each attaching point will determine the corresponding mode shapes. Since the order of [B] is 4n + 4, where n is the total number of point masses and concentrated springs, the "explicit" mathematical expression for the existing approach becomes lengthily intractable if n > 2. The "numerical assembly method"(NAM) introduced in this paper aims at improving the last drawback of the existing approach. The "exact"solutions in this paper refer to the numerical results obtained from the "continuum" models for the classical analytical approaches rather than from the "discretized" ones for the conventional finite element methods.