diff --git a/dart/api/internal/polar.py b/dart/api/internal/polar.py
index 69f5968bda90344dcc97352f842cf217c33ceb1b..bc3d3bf540360c266ada8b2af0775c9ed52bfa48 100644
--- a/dart/api/internal/polar.py
+++ b/dart/api/internal/polar.py
@@ -81,34 +81,33 @@ class Polar:
         """Compute the polar
         """
         # initialize the loop
-        ac = 0 if self.alphas[0] == self.alphas[-1] else 1
         self.Cls = []
         self.Cds = []
         self.Cms = []
         print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET)
-        for alpha in self.alphas:
+        for i in range(len(self.alphas)):
             # define current angle of attack
-            alpha = alpha*math.pi/180
-            # update problem
+            alpha = self.alphas[i]*math.pi/180
+            acs = '_{:04d}'.format(i)
+            # update problem and reset ICs to improve robustness in transonic cases
             self.pbl.update(alpha)
+            self.sol.reset()
             # run the solver and save the results
-            print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
+            print(ccolors.ANSI_BLUE + 'Running the solver for', self.alphas[i], 'degrees', ccolors.ANSI_RESET)
             self.tms["solver"].start()
             self.sol.run()
             self.tms["solver"].stop()
-            self.sol.save(self.wrtr, ac)
+            self.sol.save(self.wrtr, acs)
             # extract Cp
             if self.dim == 2:
                 Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
-                tU.write(Cp, 'Cp_airfoil_'+str(ac)+'.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
+                tU.write(Cp, f'Cp_airfoil{acs}.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
             elif self.dim == 3 and self.format == 'vtk' and self.slice:
-                dU.writeSlices(self.msh.name, self.slice, self.tag, ac)
+                dU.writeSlices(self.msh.name, self.slice, self.tag, acs)
             # extract force coefficients
             self.Cls.append(self.sol.Cl)
             self.Cds.append(self.sol.Cd)
             self.Cms.append(self.sol.Cm)
-            # end loop
-            ac+=1
 
     def disp(self):
         """Display the results and draw the polar
@@ -127,6 +126,6 @@ class Polar:
         print(self.tms)
         # plot results
         if self.alphas[0] != self.alphas[-1]:
-            tU.plot(self.alphas, self.Cls, "alpha", "Cl", "")
-            tU.plot(self.alphas, self.Cds, "alpha", "Cd", "")
-            tU.plot(self.alphas, self.Cms, "alpha", "Cm", "")
+            tU.plot(self.alphas, self.Cls, 'alpha', 'Cl', '')
+            tU.plot(self.alphas, self.Cds, 'alpha', 'Cd', '')
+            tU.plot(self.alphas, self.Cms, 'alpha', 'Cm', '')
diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py
index 74c0025882f5de8f6be7732433c8500b8a5b109d..1fb871c6daa3ec53756ba559fc13537f565706d8 100644
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
@@ -86,8 +86,9 @@ class Trim:
         while True:
             # define current angle of attack
             print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
-            # update problem
+            # update problem and reset ICs to improve robustness in transonic cases
             self.pbl.update(self.alpha)
+            self.sol.reset()
             # run the solver
             self.tms["solver"].start()
             status = self.sol.run()
diff --git a/dart/utils.py b/dart/utils.py
index 012b86c2fa0d113660640d6d6c24b5ddceada94d..4f268a98c108f79dd3d2e5ce34cac5217de997bc 100644
--- a/dart/utils.py
+++ b/dart/utils.py
@@ -111,7 +111,7 @@ def extract(vElems, vData):
         i += 1
     return data
 
-def writeSlices(mshName, ys, wId, n = 0):
+def writeSlices(mshName, ys, wId, sfx=''):
     """Write slice data for each ys coordinate along the span (only works with VTK)
 
     Parameters
@@ -122,8 +122,8 @@ def writeSlices(mshName, ys, wId, n = 0):
         y-coordinates of slices
     wId : int
         index of physical tag to slice (see gmsh)
-    n : int (optional, > 0)
-        index of mesh/solution file (will be appended to mesh file name)
+    sfx : string (optional)
+        suffix that will be appended to mesh/solution file name
     """
     import numpy as np
     try:
@@ -134,10 +134,7 @@ def writeSlices(mshName, ys, wId, n = 0):
     import tboxVtk.reader as vtkR
     import tboxVtk.cutter as vtkC
     reader = vtkR.Reader()
-    if n == 0:
-        reader.open(mshName)
-    else:
-        reader.open(mshName+'_'+str(n))
+    reader.open(mshName+sfx)
     cutter = vtkC.Cutter(reader.reader.GetOutputPort())
     for i in range(0, len(ys)):
         pts, elems, vals = cutter.extract(cutter.cut(wId, [0., ys[i], 0.], [0., 1., 0.]), 2, ['Cp'])
@@ -147,7 +144,4 @@ def writeSlices(mshName, ys, wId, n = 0):
         data = np.hstack((data, vals['Cp']))
         tU.sort(elems, data)
         data = np.vstack((data, data[0,:]))
-        if n == 0:
-            tU.write(data, 'slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
-        else:
-            tU.write(data, 'slice_'+str(n)+'_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
+        tU.write(data, mshName+sfx+'_slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')