diff --git a/pybnf/algorithms.py b/pybnf/algorithms.py
index 9f73ab97..b1300387 100644
--- a/pybnf/algorithms.py
+++ b/pybnf/algorithms.py
@@ -860,12 +860,68 @@ def run(self, log_prefix, scheduler_node=None, resume=None, debug=False):
             logger.warning('Tried to move pickled algorithm, but it was not found')
 
         if (isinstance(self, SimplexAlgorithm) or self.config.config['refine'] != 1) and self.bootstrap_number is None:
+            # Make plot if requested
+            if self.config.config['plot_results']:
+                self.plot_result(best_pset)
+
             # End of fitting; delete unneeded files
             if self.config.config['delete_old_files'] >= 1:
                 run(['rm', '-rf', self.sim_dir])
 
         logger.info("Fitting complete")
 
+    def plot_result(self, best_pset):
+        try:
+            import matplotlib.pyplot as plt
+        except ImportError:
+            logger.exception('import matplotlib.pyplot failed')
+            print1('Could not load matplotlib. Skipping plot.')
+            return
+        logger.info('Plotting the results to figure')
+        job = Job(self.model_list, best_pset, 'plot',
+                  self.sim_dir, self.config.config['wall_time_sim'], None,
+                  self.config.config['normalization'], bool(self.config.config['delete_old_files']))
+        plot_res = run_job(job, False, self.failed_logs_dir)
+
+        simdata = plot_res.simdata
+        for z, model in enumerate(simdata):
+            xlist = []
+            ylist = []
+            xpointslist = []
+            ypointslist = []
+            titlelist = []
+            xlabellist = []
+            ylabellist = []
+            for suffix in simdata[model]:
+                if suffix not in self.exp_data[model]:
+                    continue
+                sdata = simdata[model][suffix]  # A Data object
+                edata = self.exp_data[model][suffix]  # A Data object
+                for variable in sdata.cols:
+                    if variable == sdata.indvar:
+                        continue
+                    if variable in edata.cols:
+                        xlist.append(sdata[sdata.indvar])
+                        ylist.append(sdata[variable])
+                        xpointslist.append(edata[sdata.indvar])
+                        ypointslist.append(edata[variable])
+                        xlabellist.append(sdata.indvar)
+                        ylabellist.append(variable)
+                        titlelist.append(suffix)
+            cols = min(7, len(xlist))
+            rows = 1 + (len(xlist)-1) // 7
+            plt.figure(figsize=(cols*3, rows*3))
+            for i in range(len(xlist)):
+                plt.subplot(rows, cols, i+1)
+                plt.plot(xpointslist[i], ypointslist[i], 'ko', markersize=3)
+                plt.plot(xlist[i], ylist[i], 'r-')
+                plt.xlabel(xlabellist[i])
+                plt.ylabel(ylabellist[i])
+                plt.title(titlelist[i])
+            plt.suptitle(model)
+            plt.tight_layout()
+            plt.show(block=(z == len(simdata)-1))
+
     def cleanup(self):
         """
         Called before the program exits due to an exception.
diff --git a/pybnf/config.py b/pybnf/config.py
index 4cd01305..a5d43740 100644
--- a/pybnf/config.py
+++ b/pybnf/config.py
@@ -137,8 +137,7 @@ def default_config():
             'output_every': 20, 'initialization': 'lh', 'refine': 0, 'bng_command': bng_command, 'smoothing': 1,
             'backup_every': 1, 'time_course': (), 'param_scan': (), 'min_objective': -np.inf, 'bootstrap': 0,
             'bootstrap_max_obj': None, 'ind_var_rounding': 0, 'local_objective_eval': 0, 'constraint_scale': 1.0,
-            'sbml_integrator': 'cvode', 'parallel_count': None,
-
+            'sbml_integrator': 'cvode', 'parallel_count': None,'plot_results': 0,
             'mutation_rate': 0.5, 'mutation_factor': 0.5, 'islands': 1, 'migrate_every': 20, 'num_to_migrate': 3,
             'stop_tolerance': 0.002, 'de_strategy': 'rand1',
 
diff --git a/pybnf/parse.py b/pybnf/parse.py
index 68e4616e..78ee2b88 100644
--- a/pybnf/parse.py
+++ b/pybnf/parse.py
@@ -21,7 +21,7 @@
                'local_min_limit', 'reserve_size', 'burn_in', 'sample_every', 'output_hist_every',
                'hist_bins', 'refine', 'simplex_max_iterations', 'wall_time_sim', 'wall_time_gen', 'verbosity',
                'exchange_every', 'backup_every', 'bootstrap', 'crossover_number', 'ind_var_rounding',
-               'local_objective_eval', 'reps_per_beta']
+               'local_objective_eval', 'reps_per_beta', 'plot_results']
 numkeys_float = ['extra_weight', 'swap_rate', 'min_objective', 'cognitive', 'social', 'particle_weight',
                  'particle_weight_final', 'adaptive_n_max', 'adaptive_n_stop', 'adaptive_abs_tol', 'adaptive_rel_tol',
                  'mutation_rate', 'mutation_factor', 'stop_tolerance', 'step_size', 'simplex_step', 'simplex_log_step',