Light trapping in microcavities and non-homogeneous optical media has become increasingly important for the development of quantum information processing, slow-light devices, integrated photonic circuits, and bio-inspired optical systems. However, conventional optical trapping models based on fixed geometric configurations or standard nonlinear Schrödinger formulations are generally unable to accurately represent the combined effects of high-order dispersion, spatially varying material properties, and structural asymmetry encountered in realistic optical environments. To address these limitations, a hybrid shifted Legendre–classical Laguerre operational matrix method was developed for the numerical solution of the variable-coefficient fifth-order Korteweg–de Vries (KdV) equation. The spatial domain, $x\in[0,L]$, is discretized using shifted Legendre polynomials to facilitate accurate enforcement of algebraic boundary conditions, whereas the semi-infinite temporal domain, $t\in[0,\infty)$, is represented by classical Laguerre polynomials that remain orthogonal with respect to the exponentially decaying weight function $w(t)=e^{-t}$. By combining spectral collocation with Chebyshev–Gauss–Lobatto nodes, the governing nonlinear partial differential equation was transformed into a square system of nonlinear algebraic equations, which was subsequently solved using a trust-region dogleg algorithm. Numerical results show that, under the variable coefficients, the pulse is not trapped: over $t\in[0,5]$, its peak drifts from $x\approx5$ toward $x\approx8.9$ and its amplitude decays from $\approx2.0$ to $\approx0.7$ as it sheds dispersive radiation, so the monitored $L^{2}$ quantity decreases to about 0.18 of its initial value; the proposed scheme resolves this drift-and-decay evolution stably. A controlled comparison against a fourth-order finite-difference scheme under identical time integration shows that the spectral discretization reaches a given accuracy with roughly an order of magnitude fewer unknowns and is several orders of magnitude more accurate at equal degrees of freedom. The fixed-point iteration is convergent, with the Jacobian spectral radius below one. The proposed framework provides an efficient and mathematically rigorous spectral methodology for modeling variable-coefficient higher-order soliton dynamics in non-homogeneous media, and establishes a unified computational framework for such problems. Although the Laguerre basis is defined on $[0,\infty)$, all reported simulations are evaluated on the finite window $t\in[0,5]$; the semi-infinite construction supplies a spectrally accurate temporal basis with the natural weight $e^{-t}$, and $[0,5]$ covers the dynamics of interest.