We present a finite element solver for a thermodynamically consistent model of multicomponent electrolyte systems. Rooted in non-equilibrium thermodynamics, the model rigorously incorporates mass conservation, charge neutrality, and entropy production, while capturing essential physical effects such as ion solvation, finite ion size, and pressure coupling. The governing equations-comprising N-1 partial mass balances, the electrostatic Poisson equation, and a modified momentum balance-are formulated regarding electrostatic potential, atomic fractions, and pressure to ensure numerical stability and physical fidelity. The solver is implemented using the FEniCSx framework, enabling flexible simulation of one- and two-dimensional systems under complex boundary conditions. Through systematic benchmark studies and parametric analyses, we validate the model against classical Nernst-Planck formulations and highlight its advantages in avoiding non-physical concentration profiles, particularly in high-field or high-concentration regimes. Simulations explore the effects of solvation number, Debye length, compressibility, and applied electric field on double layer formation, ionic distribution, and pressure response. We publicly provide the documented and validated solver framework.