A collocation method using cubic splines is developed and applied to simulate steady and time-dependent, including turbulent, thermally convecting flows for two-dimensional compressible fluids. The state variables and the fluxes of the conserved quantities are approximated by cubic splines in both space direction. This method is shown to be nu- merically conservative and to have a local truncation error proportional to the fourth power of the grid spacing. A 'dual-staggered' Cartesian grid, where energy and momen- tum are updated on one grid and mass density on the other, is used to discretize the flux form of the compressible Navier-Stokes equations. Each grid-line is staggered so that the fluxes, in each direction, are calculated at the grid midpoints. This numerical method is validated by simulating thermally convecting flows, from steady to turbulent, reproducing known results. Once validated, the method is used to investigate many aspects of thermal convection with high numerical accuracy.