The magnetic fields generated in non-central heavy-ion collisions are among the strongest fields produced in the universe, reaching magnitudes comparable to the scale of strong interactions. Backed by model simulations, we expect the resulting field to be spatially modulated, deviating significantly from the commonly considered uniform profile. In this work, we present the next step to improve our understanding of the physics of quarks and gluons in heavy-ion collisions by adding an inhomogeneous magnetic background to our lattice QCD simulations. We simulate 2+1 staggered fermions with physical quark masses for a range of temperatures covering the QCD phase transition. We assume a $1/\cosh(x)^2$ function to model the field profile and vary its strength to analyze the impact on the chiral condensate and the Polyakov loop. These order parameters show non-trivial spatial features due to the interplay between the sea and the valence effects as the system approaches the crossover temperature. We extrapolate these quantities to the continuum limit, draw the phase diagram in the $T$-$B$ plane and interpret the implications of an inhomogeneous B to QCD physics. We also find that in this set-up, the system develops steady electric currents which flow in equilibrium. We use these currents to present our new method of obtaining the magnetic susceptibility of the QCD medium and compare it to previously established techniques.