  We have developed a new computational method to explore astrophysical and
heliophysical phenomena, especially those considerably influenced by
non-thermal energetic particles. This novel approach considers the backreaction
from these energetic particles by incorporating the non-thermal fluid pressure
into Magnetohydrodynamics (MHD) equations. The pressure of the non-thermal
fluid is evaluated from the energetic particle distribution evolved through
Parker's transport equation, which is solved using stochastic differential
equations. We implement this method in the HOW-MHD code (Seo \& Ryu 2023),
which achieves 5th-order accuracy. We find that without spatial diffusion, the
method accurately reproduces the Riemann solution in the hydrodynamic shock
tube test when including the non-thermal pressure. Solving Parker's transport
equation allows the determination of pressure terms for both relativistic and
non-relativistic non-thermal fluids with adiabatic indices
$\gamma_{\rm{NT}}=4/3$ and $\gamma_{\rm{NT}}=5/3$, respectively. The method
also successfully replicates the Magnetohydrodynamic shock tube test with
non-thermal pressure, successfully resolving the discontinuities within a few
cells. Introducing spatial diffusion of non-thermal particles leads to marginal
changes in the shock but smooths the contact discontinuity. Importantly, this
method successfully simulates the energy spectrum of the non-thermal particles
accelerated through shock, which includes feedback from the non-thermal
population. These results demonstrate that this method is very powerful for
studying particle acceleration when a significant portion of the plasma energy
is taken by energetic particles.
