This study focuses on the Rijke tube problem, which includes features relevant to the modelling of thermoacoustic coupling in reactive flows: a compact acoustic source, an empirical model for the heat source and nonlinearities. This thermoacoustic system features both linear and nonlinear flow regimes with complex dynamical behaviour. In order to synthesize accurate time series, we tackle this problem from a numerical point of view, and start by proposing a dedicated solver designed for dealing with the underlying stiffness – in particular, the retarded time and the discontinuity at the location of the heat source. Stability analysis is performed on the limit of low-amplitude disturbances by using the projection method proposed by Jarlebring (PhD thesis, Technische Universität Braunschweig, 2008), which alleviates the problems arising from linearization with respect to the retarded time. The results are then compared with the analytical solution of the undamped system and with the results obtained from Galerkin projection methods commonly used in this setting. This analysis provides insight into the consequences of the various assumptions and simplifications that justify the use of Galerkin expansions based on the eigenmodes of the unheated resonator. We demonstrate that due to the presence of a discontinuity in the spatial domain, the eigenmodes in the heated case predicted by using Galerkin expansion show spurious oscillations resulting from the Gibbs phenomenon. Finally, time series in the fully nonlinear regime, where a limit cycle is established, are analysed and dominant modes are extracted. By comparing the modes of the linear regime to those of the nonlinear regime, we are able to illustrate the mean-flow modulation and frequency switching, which appear as the nonlinearities become significant and ultimately affect the form of the limit cycle. Analysis of the saturated limit cycles shows the presence of higher-frequency modes, which are linearly stable but become significant through nonlinear growth of the signal. This bimodal effect is not exhibited when the coupling between different frequencies is not accounted for. In conclusion, a dedicated solver for capturing thermoacoustic instability is proposed and methods for analysing linear and nonlinear regions of the resulting time series are introduced.