Our main contribution is the introduction of a modified problem approach, where the interval $[0,T]$ is partitioned into n equal subintervals, and the variable-order fractional derivative $q(\ell ) $ is approximated by a piecewise constant function $q_{m} $ on these subintervals. This transforms the original variable-order initial value problem (IVP) into a sequence of fractional differential equations of constant order, for which classical existence, uniqueness, and stability results can be rigorously established. We prove that as $n \to \infty $ , the piecewise sequence $(q_{n})$ converges to the original variable order $q(\ell ) $ , and that the solutions $x_{m}(\ell ) $ converge to $x(\ell ) $ , the exact solution of the variable-order IVP. This method facilitates the direct application of classical fixed-point theorems to establish existence, uniqueness, and Ulam–Hyers stability results for nonlinear variable-order fractional differential equations, overcoming major analytical difficulties associated with variable-order operators. The approach is applied to model anomalous diffusion in heterogeneous porous media, where the variable fractional order captures the depth-dependent memory and transport properties of layered soils or biological tissues. The physical example demonstrates how the method accurately describes the transport of pollutants or chemical solutes through media with spatially varying structure, providing a robust framework for analyzing real-world diffusion phenomena in complex environments.