The two-stream instability has wide range of astrophysical applications starting from gamma-ray bursts and pulsar glitches to cosmology. We consider one dimensional weakly relativistic Zakharov equations and describe nonlinear saturation of the oscillating two-stream instability using a three dimensional dynamical system resulting form a truncation of the nonlinear Schrodinger equation to three modes. The equilibrium points of the model are determined and their stability natures are discussed. Using the tools of nonlinear dynamics such as the bifurcation diagram, Poincaré maps, and Lyapunav exponents, existence of periodic, quasi-periodic, and chaotic solutions are established in the dynamical system. Interestingly, we observe the multistable behavior in this plasma model. The system has multiple attractors depending on the initial conditions. We also notice that the relativistic parameter plays the role of control parameter in the model. The theoretical results presented in this paper may be helpful for better understanding of space and astrophysical plasmas. © 2015 AIP Publishing LLC.