Abstract Knowledge of palaeointensity variation is required for determining the full vector variation of the geomagnetic field as a function of geological time. This provides essential constraints for numerical geodynamo models. To date, most palaeointensity determination methods are laborious, characterised by rather low success rates, and demand substantial processing time. The rocks under investigation must obey stringent criteria to yield faithful palaeointensities: the magnetic particles must be single domain, the natural remanent magnetisation must be a thermoremanent magnetisation, and during the successive heating steps in the laboratory no chemical alteration should occur. Here, we describe a new method that allows all magnetic domain states to be processed, i.e. it does not require single domain particles. The method makes use of the linearity of partial thermoremanent magnetisation (pTRM) with the applied laboratory field. Multiple specimens are used so that every sample is exposed only once to a laboratory field, warranting that all samples experienced the same magnetic history. Through the limited number of thermal steps alteration effects are reduced as well. The laboratory pTRM and natural remanent magnetisation (NRM) of the specimens are oriented parallel to minimise the effects of high-temperature tails that affect multidomain minerals. The pTRM acquisition temperature is selected below the temperature at which chemical alteration sets in and above the temperature trajectory where secondary viscous NRM components occur. The procedure requires a lower number of steps than any other palaeointensity method, reducing significantly the total time needed per rock unit. We propose to name the new protocol ‘multispecimen parallel differential pTRM method’. It provides the correct answer to ∼ 5% for artificial samples and natural rocks containing multidomain magnetic particles that were given a laboratory TRM of known intensity. Application to the Paricutin September–December 1943 lava flow (three sites, 7 specimens per site) yields a weighted mean of 45.9 ± 1.25 μT, within uncertainty margin of the expected value of 45.0 μT.