It is challenging to develop static and dynamic models of structures interacting with nonlinear and viscoelastic materials; yet there are many real world examples of structures that incorporate such materials. Furthermore, it is computationally expensive to predict the steady state response of these structures to static and harmonic loads, even when using the simpler models in the literature. The example of a pinned-pinned beam interacting with polyurethane foam that reacts in both tension and compression is focused on. The steady state solution is expressed as the sum of an arbitrary number of modes and Galerkin's method is used to derive modal amplitude equations.  For polynomial-type nonlinearities, it is possible to speed up computation time by using a convolution method to evaluate integral terms in the model. Also, incremental harmonic balance is used to make the frequency response predictions more efficient. 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. Using this solution method, the influence of various parameters e.g., loading configuration, excitation amplitude, linear and nonlinear stiffness, on the response of the beam is studied.