This article presents a new form of the heat balance integral method in which the spatial variable is subdivided into equal intervals and an approximating profile is specified in each subdivision. Although the nodal refinement to the classic heat balance integral method is well documented, these studies have only been applied to classic thermal or Stefan problems with a fixed boundary condition prescribed at x = 0, and break down for more realistic boundary conditions such as constant flux, Newton cooling, or time-dependent conditions. There are also issues regarding how to start up the computation for a region which initially has zero thickness. Here we apply the so-called boundary immobilization technique, along with piecewise smooth approximating profiles, which eliminates the need for a separate starting procedure and gives far more accurate results. This is demonstrated for both thermal and Stefan problems and a variety of boundary conditions.