The Radial Basis Function - Finite Differences (RBF-FD) method presents a novel, strong-form meshless approach for solving partial differential equations. The method generalises the classical finite difference (FD) framework, enabling arbitrary domains to be discretised using homogeneously arranged collocation nodes. While the RBF-FD method performs well for problems with smooth solution fields, it may become unstable in cases involving sharp interfaces that arise from geometry or material data. Such problems commonly appear in elastoplasticity, where stress fields exhibit discontinuous differentiability across the elastic–plastic transition. To address this issue, we propose a hybrid RBF-FD formulation that combines RBF-FD with the classical FD method. The formulation utilises RBF-FD, where the stress field is discretised with augmented polyharmonic splines on virtual FD stencils prescribed to each collocation node. The divergence operator, from the balance equation, is then evaluated via the FD approach. For Neumann-type boundary-condition stabilisation, we introduce a novel stabilisation technique where evaluation points for gradient operators are shifted inside the domain. The solution procedure is implemented under a generalised plane strain assumption and verified on elastic, thermo-elasto-plastic, and visco-plastic benchmark problems. Compared to the hybrid approach, classical RBF-FD exhibits higher accuracy but reduced stability, leading to oscillatory solution fields in non-smooth elastoplastic cases. The hybrid approach effectively eliminates these oscillations while preserving the ℎ-convergence order imposed by the polynomial augmentation. In the visco-plastic cases with zero initial yield stress, both RBF-FD and the hybrid approach perform comparably due to the smooth material response. The hybrid RBF-FD approach is further applied to two complex industrial cases involving the continuous casting and cooling of steel bars on a cooling bed. A 2.5D travelling slice model is developed, demonstrating that the proposed numerical approach applies to both laboratory-scale benchmarks and also to complex industrial steel production processes.