diff --git a/pyscripts/rundt.py b/pyscripts/rundt.py
index 14e60652..4c01f9c9 100644
--- a/pyscripts/rundt.py
+++ b/pyscripts/rundt.py
@@ -16,35 +16,28 @@
 
 class UeRun():
     ''' Class containing information on run '''
-    def __init__(self, n_stor = False):
-        from time import time
+    def __init__(self, *args, n_stor = False, **kwargs):
         from numpy import array        
         from uedge import bbb, com
         # TODO: Add restore/recover from timeslice
         # TODO: Add plot timeslice directly
         # NOTE: No -> Utilize direct I/O from file instead
-        self.tstart = time()
-        self.numvar = bbb.numvar
-        self.nx = com.nx
-        self.ny = com.ny
-        self.ixpt1 = com.ixpt1[0]
-        self.ixpt2 = com.ixpt2[0]
-        self.iysptrx = com.iysptrx
         self.equationkey = array([b'te', b'ti', b'phi', b'up', b'ni', b'ng', 
             b'tg'])
 
-        self.classvars = ['slice_ni', 'slice_ng', 'slice_up', 'slice_te', 
-            'slice_ti', 'slice_tg', 'slice_phi', 'slice_dttot', 'time', 
-            'fnorm', 'nfe', 'dt_tot', 'dtreal', 'ii1', 'ii2', 'ii1fail', 
-            'ii2fail', 'dtrealfail', 'itrouble', 'troubleeq', 'troubleindex',
-            'ylfail', 'isteon', 'istion', 'isupon', 'isphion', 'isupgon',
-            'isngon', 'istgon', 'ishymol', 'nisp', 'ngsp', 'nhsp', 'nhgsp', 
-            'nzsp', 'b0', 'ncore', 'pcoree', 'pcorei', 'internaleq',
-            'internalspecies', 'yldotsfscalfail']
-
-        # Intiialize all variables to empty lists in class
-        for var in self.classvars:
-            self.__setattr__(var, [])
+        self.setupvars = {}
+        for var in [ 'numvar', 'isteon', 'istion', 'isupon', 'isphion', 
+            'isupgon', 'isngon', 'istgon', 'ishymol', 'nisp', 'ngsp', 'nhsp', 
+            'nhgsp', 'nzsp', 'b0', 'ncore', 'pcoree', 'pcorei', 'nx', 'ny',
+            'iysptrx', 'ixpt1', 'ixpt2', 'neq'
+            ]:
+            try:
+                self.setupvars[var] = getattr(bbb, var)
+            except:
+                self.setupvars[var] = getattr(com, var)
+        self.setupvars['ixpt1'] = self.setupvars['ixpt1'][0]
+        self.setupvars['ixpt2'] = self.setupvars['ixpt2'][0]        
+        super().__init__(*args, **kwargs)
 
     def itroub(self):
         ''' Function that displays information on the problematic equation '''
@@ -58,88 +51,79 @@ def itroub(self):
             'Ion momentum', 'Ion density', 'Gas density', 'Gas temperature']
 
         # Find the fortran index of the troublemaking equation
-        self.neq = bbb.neq
-        self.itrouble.append(deepcopy(argmax(abs(bbb.yldot*\
+        self.classvars['itrouble'].append(deepcopy(argmax(abs(bbb.yldot*\
             bbb.sfscal)[:bbb.neq])+1))
         print("** Fortran index of trouble making equation is:\n{}".format(\
-            self.itrouble[-1]))
+            self.classvars['itrouble'][-1]))
         # Print equation information
         print("** Number of equations solved per cell:\n    numvar = {}\n"\
-            .format(self.numvar))
+            .format(bbb.numvar))
 
-        self.troubleeq.append(mod(self.itrouble[-1]-1, bbb.numvar)+1)
+        self.classvars['troubleeq'].append(mod(self.classvars['itrouble'][-1]-1, bbb.numvar)+1)
 
         species = ''
-        self.internaleq.append([abs(x - self.itrouble[-1]).min() for x in \
+        self.classvars['internaleq'].append([abs(x - self.classvars['itrouble'][-1]).min() for x in \
             self.equations].index(0))
-        if self.equations[self.internaleq[-1]].ndim == 3:
-            self.internalspecies.append( where(\
-                self.equations[self.internaleq[-1]] == self.itrouble[-1])\
+        if self.equations[self.classvars['internaleq'][-1]].ndim == 3:
+            self.classvars['internalspecies'].append( where(\
+                self.equations[self.classvars['internaleq'][-1]] == self.classvars['itrouble'][-1])\
                 [-1][0] + 1)
-            species = ' of species {}'.format(self.internalspecies[-1])
+            species = ' of species {}'.format(self.classvars['internalspecies'][-1])
         else:
-            self.internalspecies.append(0)
+            self.classvars['internalspecies'].append(0)
 
 
         print('** Troublemaker equation is:\n{} equation{}: iv_t={}\n'\
-            .format(equationsdescription[self.internaleq[-1]], species, 
-            self.troubleeq[-1]))
+            .format(equationsdescription[self.classvars['internaleq'][-1]], species, 
+            self.classvars['troubleeq'][-1]))
 
         # Display additional information about troublemaker cell
-        self.troubleindex.append(deepcopy(bbb.igyl[self.itrouble[-1]-1,]))
-        self.dtrealfail.append(deepcopy(bbb.dtreal))
-        self.ylfail.append(deepcopy(bbb.yl[self.itrouble[-1]-1]))
-        self.yldotsfscalfail.append(deepcopy((bbb.yldot*bbb.sfscal)\
-            [self.itrouble[-1]-1]))
+        self.classvars['troubleindex'].append(deepcopy(bbb.igyl[self.classvars['itrouble'][-1]-1,]))
+        self.classvars['dtrealfail'].append(deepcopy(bbb.dtreal))
+        self.classvars['ylfail'].append(deepcopy(bbb.yl[self.classvars['itrouble'][-1]-1]))
+        self.classvars['yldotsfscalfail'].append(deepcopy((bbb.yldot*bbb.sfscal)\
+            [self.classvars['itrouble'][-1]-1]))
         print('** Troublemaker cell (ix,iy) is:\n' + \
-            '{}\n'.format(self.troubleindex[-1]))
+            '{}\n'.format(self.classvars['troubleindex'][-1]))
         print('** Timestep for troublemaker equation:\n' + \
-            '{:.4e}\n'.format(self.dtrealfail[-1]))
+            '{:.4e}\n'.format(self.classvars['dtrealfail'][-1]))
         print('** yl for troublemaker equation:\n' + \
-            '{:.4e}\n'.format(self.ylfail[-1]))
-        print('** yl*sfscal for troublemaker equation:\n' + \
-            '{:.4e}\n'.format(self.yldotsfscalfail[-1]))
+            '{:.4e}\n'.format(self.classvars['ylfail'][-1]))
+        print('** yldot*sfscal for troublemaker equation:\n' + \
+            '{:.4e}\n'.format(self.classvars['yldotsfscalfail'][-1]))
     
-    def savesuccess(self, ii1, ii2, savedir, savename, fnrm=None):
+    def savesuccess(self, savename, fnrm=None):
         from time import time
         from uedge import bbb
         from copy import deepcopy
 
-        self.time.append(time())
-        if fnrm is None:
-            bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot)
-            self.fnorm.append(deepcopy((sum((bbb.yldot[:bbb.neq]*\
-                bbb.sfscal[:bbb.neq])**2))**0.5))
-        else:
-            self.fnorm.append(fnrm)
-        self.nfe.append(deepcopy(bbb.nfe))
-        self.dt_tot.append(deepcopy(bbb.dt_tot))
-        self.dtreal.append(deepcopy(bbb.dtreal))
-        self.ii1.append(ii1)
-        self.ii2.append(ii2)
-        self.neq = bbb.neq
         try:
-            self.save('{}_UeCase.hdf5'.format(savefname.split('.')[0]))
+            self.classvars['time'].append(time())
         except:
             pass
-
-        self.save_intermediate(savedir, savename)
+        if 'fnorm' in self.classvars.keys():
+            if fnrm is None:
+                bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot)
+                self.classvars['fnorm'].append(deepcopy((sum((bbb.yldot[:bbb.neq]*\
+                    bbb.sfscal[:bbb.neq])**2))**0.5))
+            else:
+                self.classvars['fnorm'].append(fnrm)
+        for var in ['nfe', 'dt_tot', 'dtreal']:
+            if var in self.classvars:
+                self.classvars[var].append(deepcopy(getattr(bbb, var)))
+        self.save_intermediate(savename)
 
     def store_timeslice(self):
         from copy import deepcopy
         from uedge import bbb
 
-        self.slice_ni.append(deepcopy(bbb.ni))
-        self.slice_ng.append(deepcopy(bbb.ng))
-        self.slice_up.append(deepcopy(bbb.up))
-        self.slice_te.append(deepcopy(bbb.te))
-        self.slice_ti.append(deepcopy(bbb.ti))
-        self.slice_tg.append(deepcopy(bbb.tg))
-        self.slice_phi.append(deepcopy(bbb.phi))
-        self.slice_dttot.append(deepcopy(bbb.dt_tot))
+        for var in ['ni', 'ng', 'up', 'te', 'ti', 'tg', 'phi', 'dt_tot']:
+            self.classvars['slice_{}'.format(var)].append(deepcopy(getattr(bbb, var)))
 
-    def save_intermediate(self, savedir, savename):
+
+    def save_intermediate(self, savename):
         from uedge.hdf5 import hdf5_save
+        from os.path import exists
         from uedge import bbb, com
         from h5py import File
 
@@ -149,56 +133,39 @@ def save_intermediate(self, savedir, savename):
         for var in [ 'nisp', 'ngsp', 'nhsp', 'nhgsp', 'nzsp']:
             self.__setattr__(var, com.__getattribute__(var))
 
-        try:
-            hdf5_save('{}/{}_last_ii2.hdf5'.format(savedir,savename))
-        except:
-            print('Folder {} not found, saving output to cwd...'\
-                .format(savedir))
-            hdf5_save('{}_last_ii2.hdf5'.format(savename))
 
-        # Try to store ready-made Case-file, if possible
-        try:
-            self.save('{}/{}_last_ii2_Case.hdf5'.format(savedir,savename))
-        except:
-            pass
+
+        if not exists('/'.join(savename.split('/')[:-1])):
+            print('Folder {} not found, saving output to cwd...'\
+                .format(savename.split('/')[0]))
+            hdf5_save(savename)
+            
+    
 
         try:
-            file = File('{}/{}_last_ii2.hdf5'.format(savedir, savename), 'r+')
+            self.save(savename)#.replace('.hdf5', '_UeCase.hdf5'))
         except:
-            file = File('{}_last_ii2.hdf5'.format(savename), 'r+')
-
-        file.require_group('convergence')
-        group = file['convergence']
-        group.create_dataset('t_start', data=self.tstart)
-        group.create_dataset('numvar', data=self.numvar)
-        group.create_dataset('neq', data=self.neq)
-        group.create_dataset('nx', data=self.nx)
-        group.create_dataset('ny', data=self.ny)
-        group.create_dataset('ixpt1', data=self.ixpt1)
-        group.create_dataset('ixpt2', data=self.ixpt2)
-        group.create_dataset('iysptrx', data=self.iysptrx)
-        group.create_dataset('equationkey', data=self.equationkey)
-        group.create_dataset('itermx', data=self.itermx)
-        group.create_dataset('incpset', data=self.incpset)
-        group.create_dataset('ii1max', data=self.ii1max)
-        group.create_dataset('ii2max', data=self.ii1max)
-        group.create_dataset('numrevjmax', data=self.numrevjmax)
-        group.create_dataset('numfwdjmax', data=self.numfwdjmax)
-        group.create_dataset('numtotjmax', data=self.numtotjmax)
-        group.create_dataset('rdtphidtr', data=self.rdtphidtr)
-        group.create_dataset('deldt_min', data=self.deldt_min)
-        group.create_dataset('rlx', data=self.rlx)
-
-        for var in self.classvars:
-            group.create_dataset(var, data=self.__getattribute__(var))
-
-        file.close()
-    
-    def convergenceanalysis(savefname, savedir='../solutions', fig=None,
+            hdf5_save(savename)
+
+        with File(savename, 'r+') as file:
+            file.require_group('convergence')
+            group = file['convergence']
+            group.create_dataset('t_start', data=self.tstart)
+            group.create_dataset('equationkey', data=self.equationkey)
+            for var, value in self.setupvars.items():
+                group.create_dataset(var, data=value)
+            for var, value in self.classsetup.items():
+                group.create_dataset(var, data=value)
+            for var in self.classvars:
+                group.create_dataset(var, data=self.classvars[var])
+        print('Intermediate solution written to {}'.format(savename))
+
+    def convergenceanalysis(self, savefname, fig=None,
         xaxis = 'exmain', logx = False, color='k', label=None,
         ylim = (None, None)):
         from h5py import File
         from matplotlib.pyplot import subplots
+        from os.path import exists
         from datetime import timedelta
         from matplotlib.ticker import FuncFormatter
         from numpy import cumsum, ones
@@ -210,49 +177,49 @@ def convergenceanalysis(savefname, savedir='../solutions', fig=None,
                 print('Three subplots required for plots! Aborting...')
                 return
             f = fig
-        try:
-            file = File('{}/{}'.format(savedir, savefname), 'r')
-        except:
-            print('File {}/{} not found. Aborting!'.format(savedir, 
-                savefname))
+
+        if not exists(savefname):
+            print('File {} not found. Aborting!'.format(savefname))
             return
-        data = file['convergence']
-        try:
+
+        with File(savefname, 'r') as file:
             data = file['convergence']
-        except:
-            print('Convergence data not found in {}/{}. Aborting!'.format(\
-                savedir, savefname))
-            return
-        
-        if xaxis == 'exmain':
-            xlabel = 'Exmain calls'
-            xones = ones(data['ii2'][()].shape)
-            x = cumsum(xones)
-        elif xaxis == 'nfe':
-            xlabel = 'nfe internal iterations'
-            x = cumsum(data['nfe'][()][:, 0, 0])
-        elif xaxis == 'time':
-            xlabel = 'Total wall-clock time [HH:MM]'
-            x = [timedelta(t - data['t_start'][()]) for t in data['time'][()]]
-            x = data['time'][()] - data['t_start'][()]
+            try:
+                data = file['convergence']
+            except:
+                print('Convergence data not found in {}. Aborting!'.format(\
+                    savefname))
+                return
             
-        if logx is True:
-            ax[0].loglog(x, data['fnorm'][()], '-', color=color, label=label)
-            ax[1].loglog(data['dt_tot'][()], data['fnorm'][()], '-', 
-                color=color, label=label)
-            ax[2].loglog(x, data['dtreal'][()], '-', color=color, label=label)
-        else:
-            ax[0].semilogy(x, data['fnorm'][()], '-', color=color, label=label)
-            ax[1].semilogy(data['dt_tot'][()], data['fnorm'][()], '-', 
-                color=color, label=label)
-            ax[2].semilogy(x, data['dtreal'][()], '-', color=color, 
-                label=label)
+            if xaxis == 'exmain':
+                xlabel = 'Exmain calls'
+                xones = ones(data['ii2'][()].shape)
+                x = cumsum(xones)
+            elif xaxis == 'nfe':
+                xlabel = 'nfe internal iterations'
+                x = cumsum(data['nfe'][()][:, 0, 0])
+            elif xaxis == 'time':
+                xlabel = 'Total wall-clock time [HH:MM]'
+                x = [timedelta(t - data['t_start'][()]) for t in data['time'][()]]
+                x = data['time'][()] - data['t_start'][()]
+                
+            if logx is True:
+                ax[0].loglog(x, data['fnorm'][()], '-', color=color, label=label)
+                ax[1].loglog(data['dt_tot'][()], data['fnorm'][()], '-', 
+                    color=color, label=label)
+                ax[2].loglog(x, data['dtreal'][()], '-', color=color, label=label)
+            else:
+                ax[0].semilogy(x, data['fnorm'][()], '-', color=color, label=label)
+                ax[1].semilogy(data['dt_tot'][()], data['fnorm'][()], '-', 
+                    color=color, label=label)
+                ax[2].semilogy(x, data['dtreal'][()], '-', color=color, 
+                    label=label)
+            ax[1].set_title('Total exmain evaluations: {}'.format\
+                (len(data['dtreal'][()])))
 
         ax[0].set_xlabel(xlabel)
         ax[1].set_xlabel('Accumulated plasma simualtion time [s]')
         ax[2].set_xlabel(xlabel)
-        ax[1].set_title('Total exmain evaluations: {}'.format\
-            (len(data['dtreal'][()])))
         ax[0].set_ylabel('Initial fnorm')
         ax[1].set_ylabel('Initial fnorm')
         ax[2].set_ylabel('Time-step (dtreal) [s]')
@@ -269,106 +236,169 @@ def convergenceanalysis(savefname, savedir='../solutions', fig=None,
         return f
 
 
-    def failureanalysis(savefname, savedir='../solutions', equation=None):
+    def failureanalysis(self, savefname, equation=None, N=slice(None), geometry=False):
         from h5py import File
+        from os.path import exists
         from matplotlib.pyplot import subplots
         from numpy import histogram, zeros
         from matplotlib.collections import PolyCollection
         
-        f, ax = subplots(2,1, figsize=(10, 7))
-        try:
-            file = File('{}/{}'.format(savedir, savefname), 'r')
-        except:
-            print('File {}/{} not found. Aborting!'.format(savedir, 
-                savefname))
-        try:
-            data = file['convergence']
-        except:
-            print('Convergence data not found in {}/{}. Aborting!'.format(\
-                savedir, savefname))
+        # TODO: add option for plotting offending cell on geometric grid
+        if geometry is True:
+            f, ax = subplots(1,2, figsize=(12, 8), 
+                gridspec_kw={'width_ratios': [1, 1.2]})
+            f.subplots_adjust(top=0.98)
+        else:
+            f, ax = subplots(2,1, figsize=(10, 7))
+
+        if not exists(savefname):
+            print('File {} not found. Aborting!'.format(savefname))
             return
+        
+        with File(savefname, 'r') as file:
+            try:
+                data = file['convergence']
+            except:
+                print('Convergence data not found in {}. Aborting!'.format(\
+                    savefname))
+                return
 
-        if equation is not None:
-            iequation = [x.decode('UTF-8') for x in data['equationkey']]\
-                .index(equation)
-        # Bin the equation errors
-        nspecies = 1/(data['nisp'][()] + 1)
-        nbins = 7*data['nisp'][()]
-        counts, bins = histogram((data['internaleq'][()]+\
-            data['internalspecies']*nspecies)-0.5, bins=nbins, 
-            range=(-0.5,6.5))
-        h, e = histogram(data['internaleq'][()] - 0.5, bins=7, 
-            range=(-0.5,6.5))
-        ax[0].bar([x for x in range(7)], h, width=1, edgecolor='k',
-            color=(0, 87/255, 183/255))
-        ax[0].hist(bins[3*data['nisp'][()]:-1], bins[3*data['nisp'][()]:], 
-            weights=counts[3*data['nisp'][()]:], edgecolor='k', 
-            color=(255/255, 215/255, 0))
-        ax[0].set_xticks(range(7))
-        ax[0].set_xticklabels([x.decode('UTF-8') for x in \
-            data['equationkey'][()]])
-        ax[0].grid(linestyle=':', linewidth=0.5, axis='y')
-        ax[0].set_xlim((-0.5,6.5))
-        ax[0].set_ylabel('Counts')
-        for i in range(7):
-            ax[0].axvline(i-0.5, linewidth=1, color='k')
-
-        # Visualize error locations
-        nx = data['nx'][()]        
-        ny = data['ny'][()]        
-        ixpt1 = data['ixpt1'][()]
-        ixpt2 = data['ixpt2'][()]
-        iysptrx = data['iysptrx'][()]
-        frequency = zeros((nx+2, ny+2))
-
-        cells = []
-        for i in range(nx+2):
-            for j in range(ny+2):
-                cells.append([[i-.5, j-.5], [i+.5, j-.5], 
-                    [i+.5, j+.5], [i-.5, j+.5]])
-
-        polys = PolyCollection(cells, edgecolors='k', linewidth=0.5, 
-            linestyle=':')
-            
-        for i in range(len(data['itrouble'])):
-            coord = data['troubleindex'][()][i]
-            if equation is None:
-                frequency[coord[0], coord[1]] += 1
-            elif iequation == data['internaleq'][()][i]:
-                frequency[coord[0], coord[1]] += 1
-
-        polys.set_cmap('binary')
+            if equation is not None:
+                iequation = [x.decode('UTF-8') for x in data['equationkey']]\
+                    .index(equation)
+            # Bin the equation errors
+            nspecies = 1/(data['nisp'][()] + 1)
+            nbins = 7*data['nisp'][()]
+            counts, bins = histogram((data['internaleq'][N]+\
+                data['internalspecies'][N]*nspecies)-0.5, bins=nbins, 
+                range=(-0.5,6.5))
+            h, e = histogram(data['internaleq'][N] - 0.5, bins=7, 
+                range=(-0.5,6.5))
+            ax[0].bar([x for x in range(7)], h, width=1, edgecolor='k',
+                color=(0, 87/255, 183/255))
+            ax[0].hist(bins[3*data['nisp'][()]:-1], bins[3*data['nisp'][()]:], 
+                weights=counts[3*data['nisp'][()]:], edgecolor='k', 
+                color=(255/255, 215/255, 0))
+            ax[0].set_xticks(range(7))
+            ax[0].set_xticklabels([x.decode('UTF-8') for x in \
+                data['equationkey'][()]])
+            ax[0].grid(linestyle=':', linewidth=0.5, axis='y')
+            ax[0].set_xlim((-0.5,6.5))
+            ax[0].set_ylabel('Counts')
+            for i in range(7):
+                ax[0].axvline(i-0.5, linewidth=1, color='k')
+
+            # Visualize error locations
+            nx = data['nx'][()]        
+            ny = data['ny'][()]        
+            ixpt1 = data['ixpt1'][()]
+            ixpt2 = data['ixpt2'][()]
+            iysptrx = data['iysptrx'][()]
+            frequency = zeros((nx+2, ny+2))
+
+            cells = []
+            if geometry is True:
+                try:
+                    rm = file['grid/com/rm'][()] 
+                    zm = file['grid/com/zm'][()] 
+                except:
+                    raise KeyError('Grid data not found in {}'.format(\
+                        savefname))
+                for i in range(nx+2):
+                    for j in range(ny+2):
+                        cell = []
+                        for k in [1, 2, 4, 3]:
+                            cell.append((rm[i, j, k], zm[i, j, k]))
+                        cells.append(cell)
+            else:
+                for i in range(nx+2):
+                    for j in range(ny+2):
+                        cells.append([[i-.5, j-.5], [i+.5, j-.5], 
+                            [i+.5, j+.5], [i-.5, j+.5]])
+
+            polys = PolyCollection(cells, edgecolors='k', 
+                linewidth=0.5-0.45*geometry, linestyle=':')
+                
+            for i in range(len(data['itrouble'][N])):
+                coord = data['troubleindex'][N][i]
+                if equation is None:
+                    frequency[coord[0], coord[1]] += 1
+                elif iequation == data['internaleq'][N][i]:
+                    frequency[coord[0], coord[1]] += 1
+
+        if geometry is False:
+            polys.set_cmap('binary')
+        else:
+            polys.set_cmap('binary')
         polys.set_array(frequency.reshape(((nx+2)*(ny+2),)))
 
         cbar = f.colorbar(polys, ax=ax[1])
         cbar.ax.set_ylabel('N trouble'+' for {}'.format(equation)*\
             (equation is not None), va='bottom', labelpad=20)
 
-        ax[1].plot([.5, nx+.5, nx+.5, .5, .5], [.5, .5, ny+.5, ny+.5, .5], 
-            'k-', linewidth=1)
 
         ax[1].set_xlabel('Poloidal index')
         ax[1].set_ylabel('Radial index')
         ax[1].add_collection(polys)
-        ax[1].plot([.5, nx+.5],[iysptrx+.5, iysptrx+.5], 'k-', 
-            linewidth=1)
-        ax[1].plot([ixpt1+.5, ixpt1+.5], [.5, iysptrx+.5], 'k-', 
-            linewidth=1)
-        ax[1].plot([ixpt2+.5, ixpt2+.5], [.5, iysptrx+.5], 'k-', 
-            linewidth=1)
-
-        file.close()
+        if geometry is False:
+            ax[1].plot([.5, nx+.5, nx+.5, .5, .5], [.5, .5, ny+.5, ny+.5, .5], 
+                'k-', linewidth=1)
+            ax[1].plot([.5, nx+.5],[iysptrx+.5, iysptrx+.5], 'k-', 
+                linewidth=1)
+            ax[1].plot([ixpt1+.5, ixpt1+.5], [.5, iysptrx+.5], 'k-', 
+                linewidth=1)
+            ax[1].plot([ixpt2+.5, ixpt2+.5], [.5, iysptrx+.5], 'k-', 
+                linewidth=1)
+        else:
+            ax[1].set_aspect('equal')
+            ax[1].set_ylim((0.95*zm.min(), 1.05*zm.max()))
+            ax[1].set_xlim((0.95*rm.min(), 1.05*rm.max()))
+
         return f
 
 
+    def restorevalues(self):
+        ''' Restores the original UEDGE values '''
+        from uedge import bbb
+        for key, value in self.orig.items():
+            bbb.__setattr__(key, value)
+
+    def exmain_isaborted(self):
+        ''' Checks if abort is requested '''
+        from uedge import bbb
+        bbb.exmain()
+        # Abort flag set, abort case
+        if bbb.exmain_aborted == 1: 
+            # Reset flag
+            bbb.exmain_aborted == 0 
+            # Restore parameters modified by script
+            try:
+                self.restorevalues(self)
+            except:
+                self.restorevalues()
+            return True
+
+    def message(self, string, separator='-', pad='', seppad = '', 
+        nseparator=1):
+        ''' Prints formatted self.message to stdout '''
+        # TODO: add formatting for len>75 strings
+        msg = pad.strip() + ' ' + string.strip() + ' ' + pad.strip()
+        for i in range(nseparator):
+            print(seppad + separator*(len(msg)-2*len(seppad)) + seppad)
+        print(msg)
+        print(seppad + separator*(len(msg)-2*len(seppad)) + seppad)
+
+
+
 
     def converge(self, dtreal=2e-9, ii1max=5000, ii2max=5, itermx=7, ftol=1e-5,
         dt_kill=1e-14, t_stop=100, dt_max=100, ftol_min = 1e-9, incpset=7,
         n_stor=0, storedist='lin', numrevjmax=2, numfwdjmax=1, numtotjmax=0, 
         tstor=(1e-3, 4e-2), ismfnkauto=True, dtmfnk3=5e-4, mult_dt=3.4, 
-        reset=True, initjac=False, rdtphidtr=1e20, deldt_min=0.04, rlx=0.9,
+        reset=True, initjac=True, rdtphidtr=1e20, deldt_min=0.04, rlx=0.9,
 
-        tsnapshot=None, savedir='../solutions', ii2increase=1.5):
+        tsnapshot=None, savedir='../solutions', ii2increase=0, savefname=None,
+        message=None):
 
 
         ''' Converges the case by increasing dt 
@@ -440,8 +470,11 @@ def converge(self, dtreal=2e-9, ii1max=5000, ii2max=5, itermx=7, ftol=1e-5,
         from numpy import linspace, logspace, log10, append
         from copy import deepcopy
         from uedge import bbb
+        from time import time
         from os.path import exists
 
+        self.tstart = time()
+        # TODO: add kwarg savename
 
         # Check if requested save-directory exists: if not, write to cwd
         if not exists(savedir):
@@ -449,68 +482,52 @@ def converge(self, dtreal=2e-9, ii1max=5000, ii2max=5, itermx=7, ftol=1e-5,
                 savedir))
             savedir = '.'        
 
+        if savefname is None:
+            autoname = True
+            try:
+                savefname = self.casename
+            except:
+                savefname = 'dtrun'
+        else:
+            autoname = False
+        
+        savefname = '{}/{}'.format(savedir, savefname)
+
+        if autoname is True:
+            if exists('{}_last_ii2.hdf5'.format(savefname)):
+                i = 2
+                while exists('{}_{}_last_ii2.hdf5'.format(savefname, i)):
+                    i += 1
+                savefname = '{}_{}'.format(savefname, i)
+        savefname = '{}_last_ii2.hdf5'.format(savefname)
+             
+        self.classvars = {}
+        for  var in ['slice_ni', 'slice_ng', 
+            'slice_up', 'slice_te', 'slice_ti', 'slice_tg', 'slice_phi', 
+            'slice_dt_tot', 'time', 'fnorm', 'nfe', 'dt_tot', 'dtreal', 'ii1', 
+            'ii2', 'ii1fail', 'ii2fail', 'dtrealfail', 'itrouble', 'troubleeq',
+            'troubleindex', 'ylfail', 'internaleq', 'internalspecies', 
+            'yldotsfscalfail']:
+            self.classvars[var] = []
 
         self.orig = {}
-        self.orig['itermx'] = deepcopy(bbb.itermx)
-        self.orig['dtreal'] = deepcopy(bbb.dtreal)
-        self.orig['icntnunk'] = deepcopy(bbb.icntnunk)
-        self.orig['ftol'] = deepcopy(bbb.ftol)
-        self.orig['mfnksol'] = deepcopy(bbb.mfnksol)
-        self.orig['rlx'] = deepcopy(bbb.rlx)
-        self.orig['deldt'] = deepcopy(bbb.deldt)
-        self.orig['isdtsfscal'] = deepcopy(bbb.isdtsfscal)
-        self.orig['incpset'] = deepcopy(bbb.incpset)
-
+        for var in ['itermx', 'dtreal', 'icntnunk', 'ftol', 'mfnksol', 'rlx',
+            'deldt', 'isdtsfscal', 'incpset']:
+            self.orig[var] = deepcopy(getattr(bbb, var))
         if numtotjmax == 0:
             numtotjmax = numrevjmax + numfwdjmax
-
-
-        self.itermx = itermx
-        self.incpset = incpset
-        self.ii1max = ii1max
-        self.ii2max = ii2max
-        self.numrevjmax = numrevjmax
-        self.numfwdjmax = numfwdjmax
-        self.numtotjmax = numtotjmax
-        self.rdtphidtr = rdtphidtr
-        self.deldt_min = deldt_min
-        self.rlx = rlx
-
+        self.classsetup = {}
+        for var in ['itermx', 'incpset', 'ii1max', 'ii2max', 'numrevjmax',
+            'numfwdjmax', 'numtotjmax', 'rdtphidtr', 'deldt_min', 'rlx']:
+            self.classsetup[var] = locals()[var]
 # TODO: Add variable to control reduciton factor?
 # TODO: Should dtreal = min(x, t_stop) actually be t_stop or dt_max?
 
-        def restorevalues(self):
-            ''' Restores the original UEDGE values '''
-            for key, value in self.orig.items():
-                bbb.__setattr__(key, value)
-
-        def message(string, separator='-', pad='', seppad = '', 
-            nseparator=1):
-            ''' Prints formatted message to stdout '''
-            # TODO: add formatting for len>75 strings
-            message = pad.strip() + ' ' + string.strip() + ' ' + pad.strip()
-            for i in range(nseparator):
-                print(seppad + separator*(len(message)-2*len(seppad)) + seppad)
-            print(message)
-            print(seppad + separator*(len(message)-2*len(seppad)) + seppad)
-
         def scale_timestep(scaling):
             ''' Increases/decreases time-step ''' 
             bbb.dtreal *= scaling
 
-        def exmain_isaborted(self):
-            ''' Checks if abort is requested '''
-            from uedge import bbb
-            bbb.exmain()
-            # Abort flag set, abort case
-            if bbb.exmain_aborted == 1: 
-                # Reset flag
-                bbb.exmain_aborted == 0 
-                # Restore parameters modified by script
-                restorevalues(self)
-                return True
-
-        def issuccess(self, t_stop, ftol_min):
+        def issuccess(self, t_stop, ftol_min, savefname):
             ''' Checks if case is converged '''
             from datetime import timedelta
             from time import time
@@ -520,25 +537,32 @@ def issuccess(self, t_stop, ftol_min):
                 bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot)
                 self.fnrm_old = sum((bbb.yldot[:bbb.neq-1]*\
                     bbb.sfscal[:bbb.neq-1])**2)**0.5
-                self.savesuccess(ii1, ii2, savedir, bbb.label[0].strip(\
-                    ).decode('UTF-8'), self.fnrm_old)
+                self.savesuccess(savefname,
+                    #bbb.label[0].strip(\).decode('UTF-8'), 
+                    self.fnrm_old)
                 if (bbb.dt_tot>=t_stop  or  self.fnrm_old<ftol_min):
                     print('')
-                    message('SUCCESS: ' + 'fnrm < bbb.ftol'\
+                    self.message('SUCCESS: ' + 'fnrm < bbb.ftol'\
                         *(self.fnrm_old<ftol_min) + \
                         'dt_tot >= t_stop'*(bbb.dt_tot >= t_stop), pad='**', 
                         separator='*')
                     print('Total runtime: {}'.format(timedelta(
                         seconds=round(time()-self.tstart))))
-                    restorevalues(self)
+                    try:
+                        self.restorevalues(self)
+                    except:
+                        self.restorevalues()
                     return True
     
         def isfail(dt_kill):
             ''' Checks whether to abandon case '''
             if (bbb.dtreal < dt_kill):
-                message('FAILURE: time-step < dt_kill', pad='**', 
+                self.message('FAILURE: time-step < dt_kill', pad='**', 
                 separator='*')
-                restorevalues(self)
+                try:
+                    self.restorevalues(self)
+                except:
+                    self.restorevalues()
                 return True
 
         def setmfnksol(ismfnkauto, dtmfnk3):
@@ -572,22 +596,29 @@ def calc_fnrm():
         bbb.rlx = rlx
         bbb.dtreal = dtreal
         bbb.ftol = ftol
+        ii1 = 0
+        ii2 = 0
         if (bbb.iterm == 1) and (bbb.ijactot > 0):
-            message('Initial successful time-step exists', separator='')
+            self.message('Initial successful time-step exists', separator='')
         else:
-            message('Need to take initial step with Jacobian; ' + \
+            self.message('Need to take initial step with Jacobian; ' + \
                 'trying to do here', seppad='*')
             # Ensure time-step is taken
             bbb.icntnunk = 0
+            if message is not None:
+                print(message)
             # Take timestep and see if abort requested
-            if exmain_isaborted(self):
+            if self.exmain_isaborted():
                 return
             # Increase time
             # Verify time-step was successful
             if (bbb.iterm != 1):
-                restorevalues(self)
-                message('Error: converge an initial time-step first; then ' + \
-                    're-execute command', seppad='*')
+                try:
+                    self.restorevalues(self)
+                except:
+                    self.restorevalues()
+                self.message('Error: converge an initial time-step first; '
+                    'then re-execute command', seppad='*')
                 return
         bbb.incpset = incpset
         bbb.itermx = itermx
@@ -628,7 +659,7 @@ def calc_fnrm():
             bbb.dtreal = min([bbb.dtreal, dt_max]) 
             bbb.dtphi = rdtphidtr*bbb.dtreal
             bbb.deldt =  min([bbb.deldt, deldt_0, deldt_min])
-            message('Number of time-step changes = ''{} New time-step: {:.2E}\n'\
+            self.message('Number of time-step changes = ''{} New time-step: {:.2E}\n'\
                 .format((ii1+1), bbb.dtreal), pad='***', nseparator=1)
 
             # Enter for every loop except first, unless intijac == True
@@ -661,9 +692,14 @@ def calc_fnrm():
                 # Dynamically decrease ftol as the initial ftol decreases
                 bbb.ftol = max(min(ftol, 0.01*self.fnrm_old),ftol_min)
                 # Take timestep and see if abort requested
-                if exmain_isaborted(self):
+                if message is not None:
+                    print(message)
+                if self.exmain_isaborted():
                     return
-                if issuccess(self, t_stop, ftol_min):
+                if bbb.iterm == 1:
+                    self.classvars['ii1'].append(ii1)
+                    self.classvars['ii2'].append(ii2)
+                if issuccess(self, t_stop, ftol_min, savefname):
                     return
             bbb.icntnunk = 2
             bbb.isdtsfscal = 0
@@ -676,14 +712,19 @@ def calc_fnrm():
                     bbb.ftol = max(min(ftol, 0.01*self.fnrm_old),ftol_min)
                     # Take timestep and see if abort requested
 
-                    message("Inner iteration #{}".format(ii2+1), nseparator=0, 
+                    self.message("Inner iteration #{}".format(ii2+1), nseparator=0, 
 
                         separator='')
-                    if exmain_isaborted(self):
+                    if message is not None:
+                        print(message)
+                    if self.exmain_isaborted():
                         return
-                    if issuccess(self, t_stop, ftol_min):
+                    if bbb.iterm == 1:
+                        self.classvars['ii1'].append(ii1)
+                        self.classvars['ii2'].append(ii2)
+                    if issuccess(self, t_stop, ftol_min, savefname):
                         return
-                    message("Total time = {:.4E}; Timestep = {:.4E}".format(\
+                    self.message("Total time = {:.4E}; Timestep = {:.4E}".format(\
                         bbb.dt_tot-bbb.dtreal,bbb.dtreal), nseparator=0, 
                         separator='')
 # TODO: replace with more useful information
@@ -706,19 +747,528 @@ def calc_fnrm():
                 self.itroub()
                 ''' ISFAIL '''
                 if isfail(dt_kill):
-                    self.save_intermediate(savedir, bbb.label[0].strip()\
-                        .decode('UTF-8'))
+                    self.save_intermediate(savefname)
+#                        bbb.label[0].strip().decode('UTF-8'))
                     break
                 irev = 1
-                message('Converg. fails for bbb.dtreal; reduce time-step by '+\
+                self.message('Converg. fails for bbb.dtreal; reduce time-step by '+\
                     '3, try again', pad = '***', nseparator=0)
                 scale_timestep(1/(3*mult_dt))
                 bbb.dtphi = rdtphidtr*bbb.dtreal
                 bbb.deldt *=  1/(3*mult_dt) 
                 setmfnksol(ismfnkauto, dtmfnk3)
 #                bbb.iterm = 1
-#        bbb.iterm = -1 # Ensure subsequent repetitions work as intended
+        # ii1max iterations completed without convergence
+        # Indicate solver failure regardless of whether last solve was completer
+        # or not.
+        bbb.iterm = -1 # Ensure subsequent repetitions work as intended
+        # Return False to signal fail
+        return False
+
+
+    def continuation_solve(self, 
+        var, 
+        target=None, 
+        savedir=None,
+        index=None, 
+        dt=0.02,
+        dtreal=1e-5, 
+        ftol=1e-5,
+        iicont_max=7,
+        iicont_fail_max=3,
+        cutoff=1e6,
+        itermx=5,
+        incpset=8,
+        initres=1000,
+        deltafac=2,
+        dtdeltafac=2,
+        staticiter_max=3,
+        dtlim=1e4,
+        ii1max=150,
+        saveres=7,
+        **kwargs):
+        """ Solves var up to target using continuation method
+        var - string for variable, or nested dict of variables and targets
+                If nested dict, top level is variable name and subdict contains
+                keys 'target' and optionally 'index'. These values override the 
+                kwargs of the solver. If index is not set, defaults to None
+        target - target value for variable
+        index - tuple containing index of var if var is array.
+                In case of extended array slicing, e.g. setting ranges of
+                var, index must be defined as tuples of indices and/or 
+                slice objects. E.g. to set var[1:4,:,0], one would use
+                index = (slice(1,4), slice(None), 0). If ranges of var is
+                set, ensure target has the matching, appropriate dimensions
+                or is a float
+        dt - time-step used when solving using continuation
+        staticiter_max - maximum consequtive static iterations before callong 
+                time-depenent run
+        dtreal - time-step used when resorting to time-dependent solve
+        ftol - ftol to be attained
+        iicont_max - iterations in continuation loop
+        iicont_fail_max - max fails in continuation loop before going time-dependent
+        cutoff - cutoff value for required interations at current delta to attain target
+        itermx - max number of iterations to take within time-step
+        incpset - max iterations before re-calculating jacobian
+        initres - inital resolution, i.e. delta = target-var/initres
+        deltafac - factor by which delta is increased after iicont_max increases
+        dtdeltafac - factor by which delta is increased by for time-dependent run
+        dtlim - Cutoff on (target-var)/delta for which a time-dependent run is triggered
+        ii1max - maximum number of dt iterations to take
+        saveres - how many successful iterations are allowed between saves are dumped
+        kwargs passed to rundt
+        """
+        from uedge import bbb
+        from time import time
+        from os import makedirs
+        from copy import deepcopy
+        from shutil import rmtree
+        from os.path import exists
+        from numpy import array
+        from Forthon import package, packageobject
+        # TODO: icntnunk=0 on fail only?
 
+        self.tstart = time()
+        # ==== HELPERS ====
+        def changevar():
+            """ Changes var by delta """
+            setvar(min(self.lastsuccess + self.delta, 1))
+           
+        def setvar(value):
+            """ Sets var to value """
+            for key, subdict in self.var.items():
+                newvar = subdict['origvar'] + value * subdict['deltavar']
+                if subdict['index'] is None:
+                    setattr(subdict['pkgobj'], key, newvar)
+                else:
+                    getattr(subdict['pkgobj'], key)[subdict['index']] = newvar
+
+        def getvar(key, subdict):
+            """ Returns current value of var """
+            # Variable is array: only set requested index
+            if 'index' not in subdict:
+                return getattr(subdict['pkgobj'], key)
+            else:
+                try:
+                    return getattr(subdict['pkgobj'], key)[subdict['index']]
+                except:
+                    return getattr(subdict['pkgobj'], key)
+
+        def issuccess():
+            """ Checks whether case is converged at target value of var """
+            from time import time
+            from datetime import timedelta
+            self.lastsuccess += min(self.delta, 1-self.lastsuccess)
+            if self.lastsuccess >= 1:
+                print('\n===== TARGET VALUE ACHIEVED: REDUCE FNORM ====')
+                staticiter()
+                bbb.dtreal = 1e20
+
+                print('++++++++++++++++++++++++++++++++++++++++')
+                print('+++++ CONTINUATION SOLVE SUCCEEDED +++++')
+                print('++++++++++++++++++++++++++++++++++++++++')
+                print('Total runtime: {}'.format(timedelta(
+                        seconds=round(time()-self.tstart))))
+                self.savesuccess('SUCCESS_{}.hdf5'.format(self.savedir))
+                self.restorevalues()
+                return True
+            else:
+                if self.isave >= self.saveres:
+                    self.isave = 0
+                    self.savesuccess(self.savefname.format('{:.3f}'.format(self.lastsuccess).replace('.','p')))
+                else:
+                    self.isave += 1
+
+        def printinfo(printdelta=True):
+            """ Print current solve status """
+            # TODO: add printing of indices?
+            ljust = 42
+            print('    -Variable(s) being solved:')
+            for key, _ in self.var.items():
+                print(' '.ljust(ljust-3), '-',key)
+            print('{}{:.3f}%'.format('    -Progress'.ljust(ljust),
+                        self.lastsuccess*100
+            ))
+            if printdelta is not False:
+                print('{}{:.3f}%'.format('    -Advancing by'.ljust(ljust),
+                        self.delta*100
+                ))
+                print('{}{}'.format('    -Steps to target at current delta:'.ljust(ljust), 
+                    (int((1 - self.lastsuccess)/self.delta))))
+            print('')
+
+        def isdeltatoosmall():
+            """ Returns True if delta is below cutoff """
+            if self.delta < 1/self.cutoff:
+                return True
+            return False
+
+        def dostaticiter():
+            """ Returns True if delta is above dt-run threshold """
+            if self.delta > 1/self.dtlim:
+                return True
+            return False
+                    
+
+        def staticiter():
+            """ Tries to reduce initial fnrm by iterating at constant var """
+            from uedge import bbb
+            from copy import deepcopy
+            # TODO: always dump save on successful steady-state solution?
+            printinfo(False)
+            setvar(self.lastsuccess)
+            dtreal_orig = deepcopy(bbb.dtreal)
+            ftol_orig = deepcopy(bbb.ftol)
+            bbb.dtreal = dtreal_orig/100
+            while bbb.dtreal < 1e5:
+                ftol_old = deepcopy(bbb.ftol)
+                # Dynamically decreasing fnorm
+                bbb.pandf1 (-1, -1, 0, bbb.neq, 1., bbb.yl, bbb.yldot)
+                fnorm_old = (sum((bbb.yldot[:bbb.neq]*\
+                    bbb.sfscal[:bbb.neq])**2))**0.5
+                if bbb.dtreal > 1:
+                    bbb.dtreal = 1e20
+                print('\n===== STATIC ITERATION AT DTREAL={:.2e} ====='.format(bbb.dtreal))
+                bbb.ftol = max(min(bbb.ftol, 0.01*fnorm_old), 1e-9)
+                # Take a converging step
+                if self.exmain_isaborted():
+                    setvar(self.lastsuccess)
+                    bbb.dtreal = dtreal_orig
+                    bbb.ftol = ftol_orig
+                    return None
+                # Temporary fix to avoid Segfault when using icntnunk=1!
+                bbb.ijactot = 2
+                bbb.icntnunk = 1
+                bbb.ftol = ftol_old
+                bbb.dtreal = bbb.dtreal*5
+                # Assert steady-state
+                if bbb.iterm != 1:
+                    bbb.dtreal = dtreal_orig
+                    setvar(self.lastsuccess)
+                    bbb.ftol = ftol_orig
+                    print('\n===== STATIC FNRM REDUCTION FAILED =====\n')
+                    return False
+            self.savesuccess(self.savefname.format('{:.3f}_staticiter'.format(\
+                    self.lastsuccess).replace('.','p')
+                ))
+            print('===== CONVERGED AT STEADY STATE: RETURNING TO MAIN LOOP =====')
+            bbb.dtreal = dtreal_orig
+            bbb.ftol = ftol_orig
+            return True
+
+        def dtsolve(dtdeltafac):
+            """ Perform a time-dependent solve """
+            from uedge import bbb
+            printinfo(False)
+            bbb.icntnunk = 0
+            # Make backup of original values before entering dt run
+            # TODO: compress the storing/setting/resetting of vars
+            tstart_cont = deepcopy(self.tstart)
+            incpset_cont = deepcopy(bbb.incpset)
+            itermx_cont = deepcopy(bbb.itermx)
+            dtreal_cont = deepcopy(bbb.dtreal)
+            ftol_cont = deepcopy(bbb.ftol)
+            orig_cont = deepcopy(self.orig)
+            classvars_cont = deepcopy(self.classvars)
+            classsetup_cont = deepcopy(self.classsetup)
+            bbb.incpset = 5
+            bbb.itermx = 30
+            bbb.dtreal = 1e-5
+            bbb.ftol = 1e-8
+            abort = False
+            # Ensure a first time-step can be taken
+            dtdelta = self.lastsuccess + dtdeltafac/100
+            while bbb.iterm != 1:
+                dtdelta = self.lastsuccess + dtdeltafac/100
+                setvar(dtdelta)
+                if self.exmain_isaborted():
+                    setvar(self.lastsuccess)
+                    return False
+                # If the iteration does not succeed, reduce the time-step
+                if bbb.iterm != 1:
+                    # Advancing less than .5% dt is probably an indication of a failure
+                    if dtdeltafac < 0.5:
+                        abort = True
+                        break
+                    dtdeltafac /=2
+            if abort is False:
+                self.converge(dtreal=dtreal, savedir='.', 
+                    savefname=self.savefname.format('{:.3f}_dtrun'.format(\
+                    dtdelta).replace('.','p')).replace('.hdf5',''), 
+                    message='Solving for delta={:.3f}%'.format(dtdelta*100),
+                    ii1max=ii1max, **kwargs)
+                if bbb.iterm == 1:
+                    self.lastsuccess = dtdelta
+                # Ensure original values are still being kept track of
+                self.orig = orig_cont
+                bbb.incpset = incpset_cont 
+                bbb.itermx = itermx_cont 
+                bbb.dtreal = dtreal_cont 
+                bbb.ftol = ftol_cont 
+                orig_cont = deepcopy(self.orig)
+                self.classvars = classvars_cont
+                self.classsetup = classsetup_cont
+                self.tstart = tstart_cont
+                bbb.dtreal = dt
+            if bbb.iterm != 1:
+                print('=====================================')
+                print('===== CONTINUATION SOLVE FAILED =====')
+                print('=====================================')
+                setvar(self.lastsuccess)
+                print('Progress upon abortion: {:.3f}%'.format(
+                    self.lastsuccess*100)
+                )
+                return False
+            else:
+                return True
+
+        if isinstance(var, str): 
+            if target is None:
+                raise KeyError("Keyword argument 'target' must be set"
+                    " when evaluating a single variable!"
+                )
+        elif isinstance(var, dict):
+            if savedir is None: 
+                raise KeyError("Must define save directory name when "
+                    "evalauting multiple variables!"
+                )
+        self.savedir = savedir
+        if self.savedir is None:
+            self.savedir = bbb.label[0].decode('UTF-8').strip()
+            autosavedir = True
+        else:
+            autosavedir = False
+        if len(self.savedir) == 0:
+            self.savedir = '{}_solve'.format(var)
+            autosavedir = True
+        if autosavedir:
+            if exists(self.savedir):
+                i=2
+                while exists(self.savedir+'_{}'.format(i)):
+                    i += 1
+                print('Run-folder {} exists: saving to {}_{}'.format(\
+                    self.savedir, self.savedir, i)
+                )
+                self.savedir = '{}_{}'.format(self.savedir, i)
+            else:
+                print('Saving intermediate saves to {}'.format(self.savedir))
+        else:
+            if exists(self.savedir):
+                rmtree(self.savedir)
+        # Finalize savename w/ placeholder
+        self.savefname = '{}/progress{{}}.hdf5'.format(self.savedir)
+        makedirs(self.savedir)
+        # Additional variables to be stored in save files
+        self.classvars = {}
+        for dictkeys in ['delta', 'progress', 'delta_fail', 'progress_fail']:
+            self.classvars[dictkeys] = []
+        # Find package variable
+        # TODO: Use dict to set one or multiple variables instead of string?
+
+        # Set up variables that need to be accessed by helpers
+        self.dtlim = dtlim
+        self.cutoff = cutoff
+        # Set up control flags
+        start = True
+        nstaticiter = 0
+        iicont = 0
+        iicont_fail = 0
+        dtiter = 0
+        self.isave = saveres
+        self.saveres = saveres
+        self.delta = 1/initres
+        self.lastsuccess = 0
+
+        # Set up dictionary with variables and targets
+        self.var = {}
+        # Construct the variable dictionary based on the user input
+        if isinstance(var, str):
+            self.var = {    var: {
+                            'index': index,
+                            'target': target
+                    }}
+        else:
+            self.var = var
+            for key, subdict in self.var.items():
+                if 'target' not in subdict.keys():
+                    raise KeyError("Target not set for variable {}!".format(key))
+                elif 'index' not in subdict.keys():
+                    self.var[key]['index'] = None
+        # Complement the user input with gradients for each variable
+        for key, subdict in self.var.items():
+            if isinstance(subdict['target'], list):
+                subdict['target'] = array(subdict['target'])
+            for pkg in package():
+                pkgobj = packageobject(pkg)
+                if key in pkgobj.varlist():
+                    subdict['pkg'] = pkg
+                    subdict['pkgobj'] = pkgobj
+            if subdict['index'] is not None:
+                try:
+                    getattr(subdict['pkgobj'], key)[subdict['index']]
+                except:
+                    raise ValueError('{} does not accommodate requested indices!'.format(key))
+            subdict['origvar'] = deepcopy(getvar(key, subdict))
+            subdict['deltavar'] = deepcopy(subdict['target'] - subdict['origvar'])
+            try:
+                dist = abs(subdict['deltavar']).max()
+            except:
+                dist = abs(subdict['deltavar'])
+            if dist < 1e-10:
+                raise ValueError('Target equals current value for {}, aborting!'.format(key))
+                
+        self.classsetup = {}
+        for key, subdict in self.var.items():
+            self.setupvars[key] = getattr(subdict['pkgobj'], key)
+            self.classsetup['initial_{}'.format(key)] = getvar(key, subdict)
+            self.classsetup['delta_{}'.format(key)] = subdict['deltavar']
+            self.classsetup['target_{}'.format(key)] = subdict['target']
+            if isinstance(subdict['index'], (tuple, slice)):
+                self.classsetup['index_{}'.format(key)] = str(subdict['index']).encode("ascii", "ignore")
+            elif subdict['index'] is not None:
+                self.classsetup['index_{}'.format(key)] = subdict['index']
+            else:
+                self.classsetup['index_{}'.format(key)] = False
+
+        # Record original solver settings            
+        self.orig = {}
+        for ovar in ['itermx', 'dtreal', 'icntnunk', 'ftol', 'incpset',
+            'ismmaxuc', 'mmaxu'        
+        ]:
+            self.orig[ovar] = deepcopy(getattr(bbb, ovar))
+        # Take the initial time-step
+        changevar()
+        # Set remaining solver settings
+        bbb.dtreal = dt
+        bbb.ftol = ftol
+        bbb.incpset = incpset
+        bbb.itermx = itermx
+#        bbb.ismmaxuc = 0
+#        bbb.mmaxu = 70
+        if (bbb.iterm == 1) and (bbb.ijactot > 0):
+            self.message('Initial successful time-step exists', separator='')
+        else:
+            self.message('Need to take initial step with Jacobian; ' + \
+                'trying to do here', seppad='*')
+            printinfo()
+            # Ensure preconditioner is calculated
+            bbb.icntnunk = 0
+            # Take timestep and see if abort requested
+            if self.exmain_isaborted():
+                setvar(self.lastsuccess)
+                return
+            # Increase time
+            # Verify time-step was successful
+            if (bbb.iterm != 1):
+                self.restorevalues()
+                setvar(self.lastsuccess)
+                self.message('Error: converge an initial time-step first; then ' + \
+                    're-execute command', seppad='*')
+                return
+        # Start outer loop
+        while True:
+            # Ensure we don't exceed the target value
+            # Re-eval preconditioner: omit first step
+            if start is False:
+                bbb.icntnunk = 0
+            else: 
+                bbb.icntnunk = 1
+                start = False
+            # Reset counters
+            iicont = 0
+            iicont_fail = 0
+            # Start inner loop
+            while iicont < iicont_max:
+                print('===== MAIN LOOP {}/{} ====='.format(iicont+1, iicont_max))
+                printinfo()
+                # TODO: how to get intial fnrm rather than last fnorm?
+                fnorm_old = (sum((bbb.yldot[:bbb.neq]*\
+                    bbb.sfscal[:bbb.neq])**2))**0.5
+                bbb.ftol = min(max(fnorm_old/100, 1e-6), ftol)
+                # Take step and check whether abort is requested
+                if self.exmain_isaborted():
+                    setvar(self.lastsuccess)
+                    return
+                bbb.ftol = ftol
+                bbb.icntnunk = 1
+                # Failure to take time-step
+                if bbb.iterm != 1:
+                    # Check whether an excorbant number of iterations are necessary to converge: fail if yes
+                    self.classvars['delta_fail'].append(self.delta)
+                    self.classvars['progress_fail'].append(self.lastsuccess)
+                    if isdeltatoosmall():
+                        print('=====================================')
+                        print('===== CONTINUATION SOLVE FAILED =====')
+                        print('=====================================')
+                        print('Last successful step for {}: {:.4e}'.format(\
+                            self.var, self.lastsuccess)
+                        )
+                        return
+                    # Start trying to reset convergence
+                    else:
+                        # Reset success counter
+                        iicont = 0
+                        iicont_fail += 1
+                        # Delta has been reduced too many times: do some moves to knock the case loose
+                        if iicont_fail >= iicont_fail_max:
+                            print('===== RECOVERY LOOP ENTERED =====')
+                            # Iterate the case at the last successful 
+                            bbb.icntnunk = 0
+                            if dostaticiter():
+                                staticstatus = staticiter()
+                                # Flag for abortion
+                                if staticstatus is None:
+                                    return 
+                                elif staticstatus is False:
+                                    dtcall=True
+                                else:
+                                    dtcall=False
+                                nstaticiter += 1
+                            # Enter loop for time-dependent simulations
+                            if (bbb.iterm != 1) or (dtcall is True) or (nstaticiter==staticiter_max):
+                                print('===== ENTER TIME-DEPENDENT SOLVE =====')
+                                if dtsolve(dtdeltafac) is False:
+                                    return
+                                else:
+                                    # Force-break inner loop
+                                    iicont = 1e5
+                                    dtcall = False
+                                    bbb.icntnunk = 0
+                                    if issuccess():
+                                        return
+                                    changevar()
+                            # Initial fnrm successfully reduced, continue increasing var
+                            else:
+                                iicont_fail = 0
+                                changevar()
+                        # Try to attain convergence again by decreasing change in var
+                        else:
+                            print('\n===== FAIL {}/{} ====='.format(iicont_fail, iicont_fail_max))
+                            print('===== DECREASE DELTA AND TRY AGAIN =====')
+                            # Reset to last successful step
+                            setvar(self.lastsuccess)
+                            # Decrease change in variable
+                            self.delta /= deltafac
+                            # Set to new, decreased change
+                            changevar()
+                            bbb.icntnunk = 0
+                # Success: update counter and last successful value
+                else:
+                    self.classvars['delta'].append(self.delta)
+                    self.classvars['progress'].append(self.lastsuccess)
+                    if issuccess():
+                        return
+                    # Check whether to increment time-step
+                    if iicont == iicont_max-1:
+                        print('\n===== INNER LOOP COMPLETED: ADVANCING DELTA =====')
+                        nstaticiter = 0
+                        self.delta *= 1.1*deltafac
+                    else:
+                        print('\n===== SUCCESS: ADVANCING VARIABLE =====')
+                    iicont += 1
+                    changevar()
+        print('EXITED LOOP: YOU SHOULD NOT BE HERE...')
+        return
 
 def rundt(**kwargs):
     runcase=UeRun()