A nonlinear mathematical model is developed in the time domain to simulate the behaviour of two identical flexibly mounted cylinders in tandem while undergoing vortex-induced vibration (VIV). Subsequently, the model is validated and modified against experimental results. Placing an array of bluff bodies in proximity frequently happens in different engineering fields. Chimney stacks, power transmission lines and oil production risers are few engineering structures that may be impacted by VIV. The coinciding of the vibration frequency with the structure natural frequency could have destructive consequences. The main objective of this study is to provide a symplectic and reliable model capable of capturing the wake interference phenomenon. This study shows the influence of the leading cylinder on the trailing body and attempts to capture the change in added mass and damping coefficients due to the upstream wake. The model is using two coupled equations to simulate the structural response and hydrodynamic force in each of cross-flow and stream-wise directions. Thus, four equations describe the fluid–structure interaction of each cylinder. A Duffing equation describes the structural motion, and the van der Pol wake oscillator defines the hydrodynamic force. The system of equations is solved analytically. Two modification terms are added to the excitation side of the Duffing equation to adjust the hydrodynamic force and incorporate the effect of upstream wake on the trailing cylinder. Both terms are functions of upstream shedding frequency (Strouhal number). Additionally, the added mass modification coefficient is a function of structural acceleration and the damping modification coefficient is a function of velocity. The modification coefficients values are determined by curve fitting to the difference between upstream and downstream wake forces, obtained from experiments. The damping modification coefficient is determined by optimizing the model against the same set of experiments. Values of the coefficients at seven different spacings are used to define a universal function of spacing for each modification coefficient so that they can be obtained for any given distance between two cylinders. The model is capable of capturing lock-in range and maximum amplitude.