<p>We present a meshless method of lines for efficiently simulating two-dimensional nonlinear reaction–diffusion systems. Using positive-definite radial kernels, we derive approximate particular solutions via the Laplace operator and superimpose them to approximate unknown fields and their spatial derivatives without mesh generation. The resulting ordinary differential system is integrated with a high-order solver. We theoretically establish and numerically verify the method’s stability and positivity-preserving properties, ensuring dynamical consistency with the underlying PDEs. Benchmark tests on cross-diffusion Brusselator systems confirm accuracy, robustness, and geometric flexibility. Compared with leading mesh-based and meshless schemes, our approach offers a compatible and high-order framework for reaction–diffusion problems on complex domains.</p>