We study the system of interacting axions and magnetic fields in the early universe after the quantum chromodynamics phase transition, when axions acquire masses. Both axions and magnetic fields are supposed to be spatially inhomogeneous. We derive the equations for the spatial spectra of these fields, which depend on conformal time. In case of the magnetic field, we deal with the spectra of the energy density and the magnetic helicity density. The evolution equations are obtained in the closed form within the mean field approximation. We choose the parameters of the system and the initial condition which correspond to realistic primordial magnetic fields and axions. The system of equations for the spectra is solved numerically. We compare the cases of inhomogeneous and homogeneous axions. The evolution of the magnetic field in these cases is different only within small time intervals. Generally, magnetic fields are driven mainly by the magnetic diffusion. We find that the magnetic field instability takes place for the amplified initial wavefunction of the homogeneous axion. This instability is suppressed if we account for the inhomogeneity of the axion.