Buildings, railway tracks, drill strings and off-shore pipelines are all treated as structures on elastic foundations in order to study their response behavior in many engineering applications. Also, flexible polyurethane foams used for cushioning in furniture, foot-ware, and automotive industries serve as foundations, and exhibit complex nonlinear viscoelastic behavior. It is challenging to develop models of systems that include these foam-like materials and are able to predict the behaviour over a wide range of loading conditions. Even when using the simpler models commonly utilized in the literature, it is computationally expensive to predict the steady-state response of these structures to static and harmonic loads. In this work a pinned-pinned beam interacting with a viscoelastic foundation which can react both in tension and compression, or in compression alone is considered. The model developed here is capable of predicting the response to static as well as dynamic forces, whether they are concentrated or distributed. If the foundation reacts only in compression, the contact region changes with beam motion and the estimation of the unknown contact region is embedded into the iterative solution procedure. The steady-state solution is expressed as the sum of an arbitrary number of modes of an undamped pinned-pinned beam and Galerkin method is used to derive equations for the modal amplitudes. Incremental harmonic balance is used to make the steady-state frequency response predictions more efficient and a pseudo arc-length continuation technique is used to track both stable and unstable solution branches. By using these computationally efficient solution approaches, it is possible to explore a much wider variety of loading conditions and also quickly determine the number of modes required for convergence of the periodic solution. By using this solution method, the influence of various system parameters on the response of the beam is studied.