This paper presents a numerical approach to compute polar perturbations in black hole systems surrounded by generic, spherically symmetric, anisotropic matter distributions. The study focuses on gravitational wave (GW) generation and propagation from extreme mass-ratio inspirals (EMRIs) in the presence of dark matter (DM) models, including Hernquist, Navarro-Frenk-White (NFW), and Einasto profiles. The authors develop a fully numerical framework to study the effects of the surrounding environment on GW emission, which is essential for understanding the dynamics of EMRIs in astrophysical settings. The paper introduces a numerical method to solve the equations governing polar perturbations in the presence of generic DM density profiles. This method involves solving for the background spacetime metric and the associated fluid variables, and then solving a set of coupled first-order ordinary differential equations (ODEs) that describe the perturbations. The framework is validated against known results for the axial sector and extended to the polar sector. The study shows that polar fluxes are not adequately described by simple gravitational redshift effects, unlike the axial case. This indicates that polar perturbations offer a new avenue for studying black hole environments through GW observations. The results demonstrate that the energy fluxes depend on the DM halo mass and compactness, with significant differences observed for different DM profiles. The authors also present a closed-form analytic fit for the DM spike profiles derived from the adiabatic growth of a BH within the Hernquist profile. These profiles are used to compute the GW energy fluxes for various DM models and configurations. The results show that the energy fluxes vary significantly depending on the DM halo properties, with the most significant deviations from the vacuum case observed for the Einasto profile. The paper concludes that the polar fluxes cannot be adequately described by post-Newtonian (PN) expansions or simple gravitational redshift effects. This highlights the importance of a fully relativistic treatment of the environment in the study of EMRIs. The results suggest that future work should focus on developing waveform models based on these flux calculations and using them in a Bayesian framework to estimate the detectability of DM overdensities. The study also emphasizes the need for a more accurate treatment of dynamical friction and backreaction effects in the context of EMRI orbital evolution.This paper presents a numerical approach to compute polar perturbations in black hole systems surrounded by generic, spherically symmetric, anisotropic matter distributions. The study focuses on gravitational wave (GW) generation and propagation from extreme mass-ratio inspirals (EMRIs) in the presence of dark matter (DM) models, including Hernquist, Navarro-Frenk-White (NFW), and Einasto profiles. The authors develop a fully numerical framework to study the effects of the surrounding environment on GW emission, which is essential for understanding the dynamics of EMRIs in astrophysical settings. The paper introduces a numerical method to solve the equations governing polar perturbations in the presence of generic DM density profiles. This method involves solving for the background spacetime metric and the associated fluid variables, and then solving a set of coupled first-order ordinary differential equations (ODEs) that describe the perturbations. The framework is validated against known results for the axial sector and extended to the polar sector. The study shows that polar fluxes are not adequately described by simple gravitational redshift effects, unlike the axial case. This indicates that polar perturbations offer a new avenue for studying black hole environments through GW observations. The results demonstrate that the energy fluxes depend on the DM halo mass and compactness, with significant differences observed for different DM profiles. The authors also present a closed-form analytic fit for the DM spike profiles derived from the adiabatic growth of a BH within the Hernquist profile. These profiles are used to compute the GW energy fluxes for various DM models and configurations. The results show that the energy fluxes vary significantly depending on the DM halo properties, with the most significant deviations from the vacuum case observed for the Einasto profile. The paper concludes that the polar fluxes cannot be adequately described by post-Newtonian (PN) expansions or simple gravitational redshift effects. This highlights the importance of a fully relativistic treatment of the environment in the study of EMRIs. The results suggest that future work should focus on developing waveform models based on these flux calculations and using them in a Bayesian framework to estimate the detectability of DM overdensities. The study also emphasizes the need for a more accurate treatment of dynamical friction and backreaction effects in the context of EMRI orbital evolution.