We investigate the nature of subcritical, finite-amplitude, double-diffusive convection in the infinite Prandtl number regime, applicable to magma chambers and the D″-layer at the core-mantle boundary, by a two-dimensional, finite-element method based on stream-function, compositional and temperature fields. Grid refinement is used for resolving the disparately-scaled thermal and chemical boundary layers present for the large ratios of the thermal to mass diffusivity (Lewis number) characteristic of magmas. In the diffusive regime a large enoughLeis required for the establishment of steady double-diffusive convection under subcritical conditions. This criticalLevaries nonlinearly with the buoyancy ratioRp, the ratio of the chemical to thermal buoyancy. For largeLe, the steady-state heat transport in the diffusive regime depends weakly onLeand approaches that found for pure thermal convection. In accordance with steady-state boundary-layer scaling, the chemical Nusselt numberNucis found numerically to vary asNuc= 1.02Le0.49Nutfor stress-free boundaries andNuc=0.96Le0.34Nutfor rigid boundaries, withNutthe thermal Nusselt number. For larger aspect-ratios a more complicated bifurcation pattern withLeis found, with the sequence ranging from no steady states, to three steady states and then to a single elongated cell, asLeis increased. Subcritical steady-state solutions can be attained by integrating the set of time-dependent double-diffusive equations. Applications of these results to the chemical boundary layers at the core-mantle boundary would suggest the D″-layer, if it is chemically stratified there, must be a time-dependent feature. Time-dependent calculations show a strong sensitivity to the initial conditions. Subcritical convective solutions in the finger regime exhibit transitions, leading to complex time-dependent flows. The tendency to form narrow cells in the subcritical, finite-amplitude, finger regimes may account for laterally variable composition in a nearly conductive thermal state. Subcritical finger instabilities are found to be able to penetrate through the entire layer in a narrow slot, as in finite Prandtl number calculations.