Introduced is an efficient numerical method to calculate the displacements of structures involving multiple, shear-deformable beams in lengthwise nonlinear elastic contact. Applicability of the method is general, and includes structures such as stacked beams, pipes, or machine components. A particularly relevant application, which is discussed and demonstrated in this paper, is the prediction of deflections and contact conditions between the rolls and strip in rolling mills used to process flat metals. The presented nonlinear displacement calculation method involves iteration using an efficient, simplified mixed finite element formulation. Nonlinearities in the general coupled beam problem can arise either from changing contact conditions between the lengthwise-coupled shear-deformable beams (opening/closing of gaps) or from nonlinear elastic coupling-stiffness between the beams. Unlike the conventional problem of end-connected beams resting on elastic foundations, this work presents an element stiffness matrix formulation and sample nonlinear solution techniques for the problem involving multiple beams that are mutually coupled along their lengths by intervening, nonlinear elastic foundations. Important practical examples of the method are given, based on production cold rolling mill data for type 301 stainless steel on a four-high mill. The examples demonstrate efficiency and value of the nonlinear coupled beam formulation in identifying appropriate machined roll profiles to produce desired strip flatness quality. The cold rolling applications presented are characterized by nonlinear Hertz foundations with stiffness-hardening, and by discontinuous contact conditions arising from nonuniform, machined roll diameter profiles. Solution efficiencies for the new formulation are compared using modified Newton–Raphson, average stiffness, and direct substitution techniques. While the average stiffness method is relatively slower, its characteristics may make it beneficial for transfer function development in flatness control systems. Application of the method to cold rolling mills illustrates the importance of including discontinuous contact nonlinearities, and enables reduced trial and error to determine suitable roll profiles. Presented is a new efficient variational formulation for computing displacements in structures involving multiple beams laterally coupled with nonlinear elastic foundations and changing contact conditions. Practical value of the method is demonstrated for identification of machined roll profiles in cold rolling mills in order to achieve desirable flatness of rolled products.