Multistabilities were found in the ocean-atmosphere flow, in a reduced order ocean-atmosphere coupled model, when the non-linear temperature equations were solved numerically. In this paper we explain how the full non-linear Stefan-Bolzmann law was numerically implemented, and the resulting change to the system dynamics compared to the original model where these terms were linearised. Multiple stable solutions were found that display distinct oceanatmosphere flows, as well as different Lyapunov stability properties. In addition, distinct Low Frequency Variability (LFV) behaviour was observed in stable attractors. We investigated the impact on these solutions of changing the magnitude of the ocean-atmospheric coupling, as well as the atmospheric emissivity to simulate an increasing greenhouse effect. Where multistabilities exist for fixed parameters, the possibility for tipping between solutions was investigated, but tipping did not occur in this version of the model where there is a constant solar forcing. This study was undertaken using a reduced-order quasi-geostrophic