diff --git a/src/feedback/GEAR/feedback.c b/src/feedback/GEAR/feedback.c
index f4f8f30961c8eb3bebd2ef99dfeccbe183b881c7..e32c72dd4914bb5b11d73d89e39de85205966282 100644
--- a/src/feedback/GEAR/feedback.c
+++ b/src/feedback/GEAR/feedback.c
@@ -225,7 +225,13 @@ void feedback_evolve_spart(struct spart* restrict sp,
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
+
+  if (star_age_beg_step < -1e-6) {
+    error("Negative age for a star");
+  }
 #endif
+  const double star_age_beg_step_safe =
+      star_age_beg_step < 0 ? 0 : star_age_beg_step;
 
   /* Reset the feedback */
   feedback_reset_feedback(sp, feedback_props);
@@ -238,7 +244,8 @@ void feedback_evolve_spart(struct spart* restrict sp,
 
   /* Compute the stellar evolution */
   stellar_evolution_evolve_spart(sp, &feedback_props->stellar_model, cosmo, us,
-                                 phys_const, ti_begin, star_age_beg_step, dt);
+                                 phys_const, ti_begin, star_age_beg_step_safe,
+                                 dt);
 
   /* Transform the number of SN to the energy */
   sp->feedback_data.energy_ejected =
diff --git a/src/feedback/GEAR/stellar_evolution.c b/src/feedback/GEAR/stellar_evolution.c
index cae77ee33e01238200dd5bcfe60dceff9d10e3fa..2dcbc0b5ee63ce3b482b0aed25565f0ae751f066 100644
--- a/src/feedback/GEAR/stellar_evolution.c
+++ b/src/feedback/GEAR/stellar_evolution.c
@@ -187,10 +187,12 @@ void stellar_evolution_compute_discrete_feedback_properties(
     const float log_m_end_step, const float m_beg_step, const float m_end_step,
     const float m_init, const int number_snia, const int number_snii) {
 
+  /* Limit the mass within the imf limits */
+  const float m_beg_lim = min(m_beg_step, sm->imf.mass_max);
+  const float m_end_lim = max(m_end_step, sm->imf.mass_min);
+
   /* Compute the average mass */
-  const float m_avg =
-      initial_mass_function_get_integral_imf(&sm->imf, m_end_step, m_beg_step) /
-      initial_mass_function_get_integral_xi(&sm->imf, m_end_step, m_beg_step);
+  const float m_avg = 0.5 * (m_beg_lim + m_end_lim);
   const float log_m_avg = log10(m_avg);
 
   /* Compute the mass ejected */