ABSTRACT Off-centred dipole configurations have been suggested to explain different phenomena in neutron stars, such as natal kicks, irregularities in polarization of radio pulsars and properties of X-ray emission from millisecond pulsars. Here, for the first time, we model magnetothermal evolution of neutron stars with crust-confined magnetic fields and off-centred dipole moments. We find that the dipole shift decays with time if the initial configuration has no toroidal magnetic field. The decay time-scale is inversely proportional to magnetic field. The octupole moment decreases much faster than the quadrupole. Alternatively, if the initial condition includes strong dipolar toroidal magnetic field, the external poloidal magnetic field evolves from centred dipole to off-centred dipole. The surface thermal maps are very different for configurations with weak B = 1013 G and strong B = 1014 G magnetic fields. In the former case, the magnetic equator is cold while in the latter case, it is hot. We model light curves and spectra of our magnetothermal configurations. We found that in the case of cold equator, the pulsed fraction is small (below a few per cent in most cases) and spectra are well described with a single blackbody. Under the same conditions, models with stronger magnetic fields produce light curves with pulsed fraction of tens of per cent. Their spectra are significantly better described with two blackbodies. Overall, the magnetic field strength has a more significant effect on bulk thermal emission of neutron stars than does the field geometry.