diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp
index 44083788f..278dda795 100644
--- a/src/BaseBinaryStar.cpp
+++ b/src/BaseBinaryStar.cpp
@@ -2554,27 +2554,26 @@ void BaseBinaryStar::EmitGravitationalWave(const double p_Dt) {
 /* 
  * Choose a timestep based on the parameters of the binary.
  *
- * This function will return the minimum of (i) a timestep based on the
- * orbital timescale of the binary and (ii) (if configured to emit GWs)
- * a timestep based on the magnitude of gravitational radiation.
+ * Returns a timestep based on the orbital timescale of the binary or, if emitting GWs,
+ * magnitude of gravitational radiation
  * 
  *
- * double ChooseTimestep(const double p_Dt)
+ * double ChooseTimestep(const double p_Multiplier)
  * 
- * @param   [IN]    p_Dt                        previous timestep in Myr
- * @return                                      new timestep in Myr
+ * @param   [IN]    p_Multiplier                timestep multiplier
+ * @return                                      timestep in Myr
  */
-double BaseBinaryStar::ChooseTimestep(const double p_Dt) {
+double BaseBinaryStar::ChooseTimestep(const double p_Multiplier) {
 
-    double newDt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep());        			// new timestep
+    double dt = std::min(m_Star1->CalculateTimestep(), m_Star2->CalculateTimestep());                   // timestep based on orbital timescale
 
-    if (OPTIONS->EmitGravitationalRadiation()) {                                                                        // emitting GWs?
-        newDt = std::min(newDt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW);                                                  // reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission
+    if (OPTIONS->EmitGravitationalRadiation()) {                                                        // emitting GWs?
+        dt = std::min(dt, -1.0E-2 * m_SemiMajorAxis / m_DaDtGW);                                        // yes - reduce timestep if necessary to ensure that the orbital separation does not change by more than ~1% per timestep due to GW emission
     }
 
-    newDt *= OPTIONS->TimestepMultiplier();	
+    dt *= p_Multiplier;	
 
-    return std::max(std::round(newDt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);                // quantised and not less than minimum
+    return std::max(std::round(dt / TIMESTEP_QUANTUM) * TIMESTEP_QUANTUM, NUCLEAR_MINIMUM_TIMESTEP);    // quantised and not less than minimum
 }
 
 
@@ -2780,7 +2779,28 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
             // evolve the current binary up to the maximum evolution time (and number of steps)
 
             double dt;                                                                                                                  // timestep
-            unsigned long int stepNum = 0;                                                                                              // initialise step number
+            if (usingProvidedTimesteps) {                                                                                               // user-provided timesteps?
+                // get new timestep
+                //   - don't quantise
+                //   - don't apply timestep multiplier
+                // (we assume user wants the timesteps in the file)
+                // 
+                // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP?
+                dt = timesteps[0];
+            }
+            else {                                                                                                                      // no - not using user-provided timesteps
+                // if user selects to emit GWs, calculate the effects of radiation
+                //     - note that this is placed before the call to ChooseTimestep() because when
+                //       emitting GWs the timestep is a function of graviational radiation
+                if (OPTIONS->EmitGravitationalRadiation()) {
+                    CalculateGravitationalRadiation();
+                }
+
+                // we want the first timestep to be small - calculate timestep and divide by 1000.0
+                dt = ChooseTimestep(OPTIONS->TimestepMultiplier() / 1000.0);                                                            // calculate timestep - make first step small
+            }
+
+            unsigned long int stepNum = 1; 
 
             while (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                     // perform binary evolution - iterate over timesteps until told to stop
 
@@ -2819,6 +2839,10 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                 }
                 else {                                                                                                                  // continue evolution
 
+                    if (OPTIONS->EmitGravitationalRadiation()) {                                                                        // emitting GWs?
+                        EmitGravitationalWave(dt);                                                                                      // yes - emit graviataional wave
+                    }
+
                     (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::POST_STELLAR_TIMESTEP);                                   // print (log) detailed output
 
                     if (OPTIONS->RLOFPrinting()) StashRLOFProperties(MT_TIMING::PRE_MT);                                                // stash properties immediately pre-Mass Transfer 
@@ -2900,6 +2924,33 @@ EVOLUTION_STATUS BaseBinaryStar::Evolve() {
                     SHOW_WARN(ERROR::TIMESTEPS_EXHAUSTED);                                                                              // show warning
                 }
 
+                if (evolutionStatus == EVOLUTION_STATUS::CONTINUE) {                                                                    // continue evolution?
+
+                    // if user selects to emit GWs, calculate the effects of radiation
+                    //     - note that this is placed before the call to ChooseTimestep() because when
+                    //       emitting GWs the timestep is a function of graviational radiation                    
+                    if (OPTIONS->EmitGravitationalRadiation()) {
+                        CalculateGravitationalRadiation();
+                    }
+
+                    m_Star2->UpdatePreviousTimestepDuration();
+                    m_Star1->UpdatePreviousTimestepDuration();
+                
+                    if (usingProvidedTimesteps) {                                                                                       // user-provided timesteps?
+                        // select a timestep
+                        //   - don't quantise
+                        //   - don't apply timestep multiplier
+                        // (we assume user wants the timesteps in the file)
+                        // 
+                        // Open question: should we clamp this to NUCLEAR_MINIMUM_TIMESTEP?
+                        dt = timesteps[stepNum];
+                    }
+                    else {                                                                                                              // no - not using user-provided timesteps
+                        dt = ChooseTimestep(OPTIONS->TimestepMultiplier());
+                    }
+
+                    stepNum++;                                                                                                          // increment stepNum
+                }
             }
 
             if (usingProvidedTimesteps && timesteps.size() > stepNum) {                                                                 // all user-defined timesteps consumed?
diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h
index 12c2fc5b8..00589179e 100644
--- a/src/BaseBinaryStar.h
+++ b/src/BaseBinaryStar.h
@@ -419,7 +419,7 @@ class BaseBinaryStar {
     void    CalculateGravitationalRadiation();
     void    EmitGravitationalWave(const double p_Dt);
 
-    double  ChooseTimestep(const double p_Dt);
+    double  ChooseTimestep(const double p_Multiplier);
 
     void    CalculateEnergyAndAngularMomentum();