This study investigates the \((2+1)\) -dimensional integro-differential Jaulent–Miodek (JM) evolution equation, a fundamental model associated with energy-dependent Schrödinger potentials and widely applicable in fluid dynamics, optics, condensed matter physics, and engineering systems. By employing the modified extended direct algebraic Method, we derive a comprehensive family of exact optical soliton solutions expressed through hyperbolic, trigonometric, exponential, Jacobi elliptic, Weierstrass elliptic, and rational functions. The obtained results reveal diverse nonlinear structures, including dark and mixed-kink solitons, breathers, periodic and singular kink solitons. These solutions are presented in parametric form and systematically visualized through 2D, 3D with contour, and density plots to elucidate their qualitative dynamics. To further examine stability and robustness, the Levenberg–Marquardt artificial neural network (LM-ANN) framework is employed, achieving high accuracy in reproducing the analytical profiles with convergence validated by fitness and regression metrics. Additionally, the numerical stability of the obtained soliton families is further verified using the Vakhitov–Kolokolov criterion, which confirms the persistence of the localized wave structures under small perturbations. The proposed hybrid framework not only advances the theoretical understanding of nonlinear wave propagation in higher-dimensional systems but also bridges classical soliton theory with modern machine learning methodologies. These findings offer promising implications for optical communication, ultrafast photonics, and nonlinear signal processing.