Determination of a magnetic field structure on a neutron star (NS) surface is an important problem of a modern astrophysics. In a presence of strong magnetic fields a thermal conductivity of a degenerate matter is anisotropic. In this paper we present 3D anisotropic heat transfer simulations in outer layers of magnetized NSs, and construct synthetic thermal light curves. We have used a different from previous works tensorial thermal conductivity coefficient of electrons, derived from the analytical solution of the Boltzmann equation by the Chapman-Enskog method. We have obtained a NS surface temperature distribution in presence of dipole-plus-quadrupole magnetic fields. We consider a case, in which magnetic axes of a dipole and quadrupole components of the magnetic field are not aligned. To examine observational manifestations of such fields we have generated thermal light curves for the obtained temperature distributions using a composite black-body model. It is shown, that the simplest (only zero-order spherical function in quadrupole component) non-coaxial dipole-plus-quadrupole magnetic field distribution can significantly affect the thermal light curves, making pulse profiles non-symmetric and amplifying pulsations in comparison to the pure-dipolar field.