Dynamic response of structures can be determined by mode superposition method which uses mode shape functions and modal coordinate functions. If the structure can be modelled with a uniform Euler–Bernoulli beam and a combination of clamped, free, pinned or sliding boundary conditions, analytical expressions for mode shape functions and modal coordinate functions are used. If analytical expressions for higher-order mode shapes or analytical expressions for the modal coordinate functions at high value of time are numerically evaluated, errors occur since the programming language cannot recognise the resulting value as a number, due to ill-conditioned analytical expressions. An analytical expression is ill-conditioned if small errors in the data may produce large errors in the solution. In this paper, we present a procedure for the exact numerical evaluation of analytical expressions for modal coordinates at high value of time while existing modified analytical expressions are used for exact calculation of higher-order mode shapes of the Euler–Bernoulli beam. The proposed procedure can be applied in the design of structure elements whose response is calculated with mode superposition method, e.g. for designing of metal frames for sawing machines in woodworking industry.