diff --git a/BayesianCalibration_spm/makeData.py b/BayesianCalibration_spm/makeData.py
index 964c151..94c74ff 100644
--- a/BayesianCalibration_spm/makeData.py
+++ b/BayesianCalibration_spm/makeData.py
@@ -41,7 +41,6 @@
 # from forwardPass import from_param_list_to_dict, from_param_list_to_str
 from forwardPass import from_param_list_to_str
 from plotsUtil_batt import *
-from prettyPlot.plotsUtil import pretty_labels
 
 # Read command line arguments
 args_spm = argument.initArg()
diff --git a/pinn_spm_param/integration_spm/main.py b/pinn_spm_param/integration_spm/main.py
index 026a504..fd8d09c 100644
--- a/pinn_spm_param/integration_spm/main.py
+++ b/pinn_spm_param/integration_spm/main.py
@@ -4,7 +4,6 @@
 import argument
 import numpy as np
 from plotsUtil_batt import *
-from prettyPlot.plotsUtil import pretty_labels, pretty_legend
 from thermo import grad_ds_a_cs_a, grad_ds_c_cs_c
 
 # Read command line arguments
@@ -135,25 +134,28 @@
 
 
 # Plot cs
-plotData(
-    [config["r_a"], config["r_c"]],
-    [sol["cs_a"] / params["csanmax"], sol["cs_c"] / params["cscamax"]],
-    params["tmax"],
-    [r"", r""],
-    [r"$c_{s,an}/c_{s,an,max}$", r"$c_{s,ca}/c_{s,ca,max}$"],
-    ["r [m]", "r [m]"],
+pretty_multi_contour(
+    listDatax=[config["r_a"], config["r_c"]],
+    listData=[
+        sol["cs_a"] / params["csanmax"],
+        sol["cs_c"] / params["cscamax"],
+    ],
+    ybound=[0, params["tmax"]],
+    listCBLabel=[r"", r""],
+    listTitle=[r"$c_{s,an}/c_{s,an,max}$", r"$c_{s,ca}/c_{s,ca,max}$"],
+    listXAxisName=["r [m]", "r [m]"],
     vminList=[0, 0],
     vmaxList=[1, 1],
 )
 
 # Plot cs
-plotData(
-    [config["r_a"], config["r_c"]],
-    [sol["cs_a"], sol["cs_c"]],
-    params["tmax"],
-    [r"", r""],
-    [r"$c_{s,an}$", r"$c_{s,ca}$"],
-    ["r [m]", "r [m]"],
+pretty_multi_contour(
+    listDatax=[config["r_a"], config["r_c"]],
+    listData=[sol["cs_a"], sol["cs_c"]],
+    ybound=[0, params["tmax"]],
+    listCBLabel=[r"", r""],
+    listTitle=[r"$c_{s,an}$", r"$c_{s,ca}$"],
+    listXAxisName=["r [m]", "r [m]"],
 )
 
 if args.verbose:
diff --git a/pinn_spm_param/integration_spm/spm_int.py b/pinn_spm_param/integration_spm/spm_int.py
index 0bb4d5b..43e78ee 100644
--- a/pinn_spm_param/integration_spm/spm_int.py
+++ b/pinn_spm_param/integration_spm/spm_int.py
@@ -1,11 +1,8 @@
 import sys
-
-import numpy as np
-
-sys.path.append("../util")
 import time
 
-from myProgressBar import printProgressBar
+import numpy as np
+from prettyPlot.progressBar import print_progress_bar
 
 
 def get_r_domain(n_r, params):
@@ -188,7 +185,7 @@ def integration(
     dR_c = config["dR_c"]
 
     if verbose:
-        printProgressBar(
+        print_progress_bar(
             0,
             n_t,
             prefix=f"Step= {0} / {n_t}",
@@ -429,7 +426,7 @@ def integration(
             )
 
         if verbose:
-            printProgressBar(
+            print_progress_bar(
                 i_t,
                 n_t - 1,
                 prefix=f"Step= {i_t} / {n_t}",
diff --git a/pinn_spm_param/main.py b/pinn_spm_param/main.py
index 1dfcd1e..1827ab0 100644
--- a/pinn_spm_param/main.py
+++ b/pinn_spm_param/main.py
@@ -18,8 +18,6 @@
     initialize_params_from_inpt,
 )
 from myNN import *
-from myparser import parseInputFile
-from myProgressBar import printProgressBar
 
 
 def do_training_only(input_params, nn):
diff --git a/pinn_spm_param/postProcess/computeManyMetrics.py b/pinn_spm_param/postProcess/computeManyMetrics.py
index 08d92ef..487f214 100644
--- a/pinn_spm_param/postProcess/computeManyMetrics.py
+++ b/pinn_spm_param/postProcess/computeManyMetrics.py
@@ -9,7 +9,7 @@
 
 import argument
 from computeError import *
-from prettyPlot.plotsUtil import plt, pretty_labels, pretty_legend
+from plotsUtil_batt import *
 from scipy import stats
 
 print("\n\nINFO: COMPUTING MANY METRICS\n\n")
diff --git a/pinn_spm_param/postProcess/plotCorrelationPINNvsData.py b/pinn_spm_param/postProcess/plotCorrelationPINNvsData.py
index b4044f8..219f65d 100644
--- a/pinn_spm_param/postProcess/plotCorrelationPINNvsData.py
+++ b/pinn_spm_param/postProcess/plotCorrelationPINNvsData.py
@@ -8,7 +8,7 @@
 import argument
 import tensorflow as tf
 from myNN import *
-from prettyPlot.plotsUtil import plt, pretty_labels, pretty_legend
+from plotsUtil_batt import *
 from tensorflow import keras
 from tensorflow.keras import layers, regularizers
 
diff --git a/pinn_spm_param/postProcess/plotData.py b/pinn_spm_param/postProcess/plotData.py
index 9e11d6c..97f3cf2 100644
--- a/pinn_spm_param/postProcess/plotData.py
+++ b/pinn_spm_param/postProcess/plotData.py
@@ -82,13 +82,13 @@ def plot_pde_data(args):
         plt.savefig(os.path.join(figureFolder, "PhiData.png"))
 
     # Plot cs
-    plotData(
-        [r_a, r_c],
-        [cs_a, cs_c],
-        params["tmax"],
-        [r"$[kmol/m^3]$", r"$[kmol/m^3]$"],
-        [r"$c_{s,an}$", r"$c_{s,ca}$"],
-        ["r [m]", "r [m]"],
+    pretty_multi_contour(
+        listDatax=[r_a, r_c],
+        listData=[cs_a, cs_c],
+        ybound=[0, params["tmax"]],
+        listCBLabel=[r"$[kmol/m^3]$", r"$[kmol/m^3]$"],
+        listTitle=[r"$c_{s,an}$", r"$c_{s,ca}$"],
+        listXAxisName=["r [m]", "r [m]"],
     )
     if not args.verbose:
         plt.savefig(os.path.join(figureFolder, "csData.png"))
diff --git a/pinn_spm_param/postProcess/plotManyLosses.py b/pinn_spm_param/postProcess/plotManyLosses.py
index 9a855fe..6290279 100644
--- a/pinn_spm_param/postProcess/plotManyLosses.py
+++ b/pinn_spm_param/postProcess/plotManyLosses.py
@@ -5,7 +5,7 @@
 
 sys.path.append("../util")
 import argument
-from prettyPlot.plotsUtil import plt, pretty_labels, pretty_legend
+from plotsUtil_batt import *
 
 print("\n\nINFO: PLOTTING MANY LOSSES\n\n")
 
diff --git a/pinn_spm_param/postProcess/plotPINNResult.py b/pinn_spm_param/postProcess/plotPINNResult.py
index 091d026..edfb4c1 100644
--- a/pinn_spm_param/postProcess/plotPINNResult.py
+++ b/pinn_spm_param/postProcess/plotPINNResult.py
@@ -154,13 +154,13 @@ def plot_pinn_result(args):
         field_dat=field_dat,
         verbose=args.verbose,
     )
-    plotData(
-        [r_test_a[:, :, 0, 0], r_test_c[:, :, 0, 0]],
-        [cs_a[:, :, 0, 0], cs_c[:, :, 0, 0]],
-        tmax,
-        [r"$[kmol/m^3]$", r"$[kmol/m^3]$"],
-        [r"$c_{s,an}$", r"$c_{s,ca}$"],
-        ["r [m]", "r [m]"],
+    pretty_multi_contour(
+        listDatax=[r_test_a[:, :, 0, 0], r_test_c[:, :, 0, 0]],
+        listData=[cs_a[:, :, 0, 0], cs_c[:, :, 0, 0]],
+        ybound=[0, tmax],
+        listCBLabel=[r"$[kmol/m^3]$", r"$[kmol/m^3]$"],
+        listTitle=[r"$c_{s,an}$", r"$c_{s,ca}$"],
+        listXAxisName=["r [m]", "r [m]"],
     )
     if not args.verbose:
         plt.savefig(
diff --git a/pinn_spm_param/postProcess/plotPINNResult_movie.py b/pinn_spm_param/postProcess/plotPINNResult_movie.py
index 8af2ea3..7a2a9f2 100644
--- a/pinn_spm_param/postProcess/plotPINNResult_movie.py
+++ b/pinn_spm_param/postProcess/plotPINNResult_movie.py
@@ -125,13 +125,13 @@ def makePlot(movieDir, solDict, params_list, indMov, stepID, data_dict=None):
         field_dat=field_dat,
         verbose=False,
     )
-    plotData(
-        [r_test_a[:, :, 0, 0], r_test_c[:, :, 0, 0]],
-        [cs_a[:, :, 0, 0], cs_c[:, :, 0, 0]],
-        tmax,
-        [r"$[kmol/m^3]$", r"$[kmol/m^3]$"],
-        [r"$c_{s,an}$", r"$c_{s,ca}$"],
-        ["r [m]", "r [m]"],
+    pretty_multi_contour(
+        listDatax=[r_test_a[:, :, 0, 0], r_test_c[:, :, 0, 0]],
+        listData=[cs_a[:, :, 0, 0], cs_c[:, :, 0, 0]],
+        ybound=[0, tmax],
+        listCBLabel=[r"$[kmol/m^3]$", r"$[kmol/m^3]$"],
+        listTitle=[r"$c_{s,an}$", r"$c_{s,ca}$"],
+        listXAxisName=["r [m]", "r [m]"],
     )
     plt.savefig(
         os.path.join(figureFolder, f"cs2D{string_params}_{indMov}.png")
diff --git a/pinn_spm_param/postProcess/plotResidualVariation.py b/pinn_spm_param/postProcess/plotResidualVariation.py
index 5f10265..f6a4dff 100644
--- a/pinn_spm_param/postProcess/plotResidualVariation.py
+++ b/pinn_spm_param/postProcess/plotResidualVariation.py
@@ -6,7 +6,7 @@
 
 sys.path.append("../util")
 import argument
-from prettyPlot.plotsUtil import plt, pretty_labels, pretty_legend
+from plotsUtil_batt import *
 
 print("\n\nINFO: PLOTTING LOSSES\n\n")
 
diff --git a/pinn_spm_param/util/generateOCP_poly_mon.py b/pinn_spm_param/util/generateOCP_poly_mon.py
index 1a013de..6d76dd9 100644
--- a/pinn_spm_param/util/generateOCP_poly_mon.py
+++ b/pinn_spm_param/util/generateOCP_poly_mon.py
@@ -5,7 +5,7 @@
 import numpy as np
 import pandas as pd
 import tf_lineInterp as tfli
-from myProgressBar import printProgressBar
+from prettyPlot.progressBar import print_progress_bar
 from scipy import odr
 from scipy.optimize import minimize
 
@@ -47,7 +47,7 @@ def polyOpt(order, xd, yd):
 
 def iterativeOpt(order_list, xd, yd):
     error = []
-    printProgressBar(
+    print_progress_bar(
         0,
         len(order_list) - 1,
         prefix=f"Poly Order = {0} / {len(order_list)-1}",
@@ -59,7 +59,7 @@ def iterativeOpt(order_list, xd, yd):
         pars = res.x
         poly = np.poly1d(pars)
         error.append(np.sum(abs(poly(xd) - yd)))
-        printProgressBar(
+        print_progress_bar(
             order,
             len(order_list) - 1,
             prefix=f"Poly Order = {order} / {len(order_list)-1}",
diff --git a/pinn_spm_param/util/init_pinn.py b/pinn_spm_param/util/init_pinn.py
index 5eab9a1..2aea885 100644
--- a/pinn_spm_param/util/init_pinn.py
+++ b/pinn_spm_param/util/init_pinn.py
@@ -13,8 +13,7 @@
 os.environ["TF_CPP_MIN_LOG_LEVEL"] = "2"
 import tensorflow as tf
 from myNN import *
-from myparser import parseInputFile
-from myProgressBar import printProgressBar
+from prettyPlot.parser import parse_input_file
 
 
 def absolute_path_check(path):
@@ -289,7 +288,7 @@ def initialize_params_from_inpt(inpt):
 
 
 def initialize_params(args):
-    inpt = parseInputFile(args.input_file)
+    inpt = parse_input_file(args.input_file)
     return initialize_params_from_inpt(inpt)
 
 
diff --git a/pinn_spm_param/util/myNN.py b/pinn_spm_param/util/myNN.py
index 3cc917d..4665a9c 100644
--- a/pinn_spm_param/util/myNN.py
+++ b/pinn_spm_param/util/myNN.py
@@ -18,7 +18,7 @@
 from custom_activations import swish_activation
 from dataTools import checkDataShape, completeDataset
 from eager_lbfgs import Struct, lbfgs
-from myProgressBar import printProgressBar
+from prettyPlot.progressBar import print_progress_bar
 from tensorflow.keras import backend as K
 from tensorflow.keras import (
     initializers,
@@ -2313,7 +2313,7 @@ def train(
         print("Using collocation points: " + self.collocationMode)
 
         # Train
-        printProgressBar(
+        print_progress_bar(
             0,
             self.nEpochs,
             prefix="Loss=%s  Epoch= %d / %d " % ("?", 0, self.nEpochs),
@@ -2776,7 +2776,7 @@ def train(
                     ) / (step + 1)
 
                     if not optimized and step % self.freq == 0:
-                        printProgressBar(
+                        print_progress_bar(
                             step,
                             self.n_batch,
                             prefix="Loss=%.2f i=%.2f b=%.2f d=%.2f r=%.2f, t/step=%.2g ms,  Epoch= %d / %d "
@@ -2828,7 +2828,7 @@ def train(
                             self.reg_col_weights,
                         )
 
-                printProgressBar(
+                print_progress_bar(
                     epoch + 1,
                     self.nEpochs,
                     prefix="Loss=%.2f i=%.2f b=%.2f d=%.2f r=%.2f, t/step=%.2g ms,  Epoch= %d / %d "
diff --git a/pinn_spm_param/util/myProgressBar.py b/pinn_spm_param/util/myProgressBar.py
deleted file mode 100644
index eb8d050..0000000
--- a/pinn_spm_param/util/myProgressBar.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import sys
-
-
-# Adapted from https://stackoverflow.com/questions/3173320/text-progress-bar-in-terminal-with-block-characters/34325723#34325723
-# Print iterations progress
-def printProgressBar(
-    iteration,
-    total,
-    prefix="",
-    suffix="",
-    decimals=1,
-    length=100,
-    fill="█",
-    printEnd="\r",
-):
-    """
-    Call in a loop to create terminal progress bar
-    @params:
-        iteration   - Required  : current iteration (Int)
-        total       - Required  : total iterations (Int)
-        prefix      - Optional  : prefix string (Str)
-        suffix      - Optional  : suffix string (Str)
-        decimals    - Optional  : positive number of decimals in percent complete (Int)
-        length      - Optional  : character length of bar (Int)
-        fill        - Optional  : bar fill character (Str)
-        printEnd    - Optional  : end character (e.g. "\r", "\r\n") (Str)
-    """
-    if total < 1e-12:
-        return
-    percent = ("{0:." + str(decimals) + "f}").format(
-        100 * (iteration / float(total))
-    )
-    filledLength = int(length * iteration // total)
-    bar = fill * filledLength + "-" * (length - filledLength)
-    print(f"\x1b[2K\r{prefix} |{bar}| {percent}% {suffix}", end=printEnd)
-    sys.stdout.flush()
-    # Print New Line on Complete
-    if iteration == total:
-        print()
-        sys.stdout.flush()
diff --git a/pinn_spm_param/util/myparser.py b/pinn_spm_param/util/myparser.py
deleted file mode 100644
index 19532fb..0000000
--- a/pinn_spm_param/util/myparser.py
+++ /dev/null
@@ -1,20 +0,0 @@
-import os
-import sys
-
-
-def parseInputFile(input_filename):
-    if not os.path.isfile(input_filename):
-        print("ERROR: No input file name found, assuming it is 'input'")
-        sys.exit()
-
-    # ~~~~ Parse input
-    inpt = {}
-    f = open(input_filename)
-    data = f.readlines()
-    for line in data:
-        if ":" in line:
-            key, value = line.split(":")
-            inpt[key.strip()] = value.strip()
-    f.close()
-
-    return inpt
diff --git a/pinn_spm_param/util/plotsUtil_batt.py b/pinn_spm_param/util/plotsUtil_batt.py
index 1e34a0c..70c64b3 100644
--- a/pinn_spm_param/util/plotsUtil_batt.py
+++ b/pinn_spm_param/util/plotsUtil_batt.py
@@ -6,6 +6,7 @@
     plt,
     pretty_labels,
     pretty_legend,
+    pretty_multi_contour,
 )
 
 
@@ -317,152 +318,6 @@ def plotField_single(
     return
 
 
-def plotData(
-    listDatax,
-    listData,
-    tmax,
-    listCBLabel,
-    listTitle,
-    listXAxisName=None,
-    vminList=None,
-    vmaxList=None,
-    globalTitle=None,
-):
-    lim = -1
-    lim_vmax_t = -1
-    lim_vmax_x = -1
-    lim_plot = -1
-    fig, axs = plt.subplots(1, len(listData), figsize=(len(listData) * 3, 4))
-    if len(listData) == 1:
-        i_dat = 0
-        data = listData[i_dat]
-        data_x = np.squeeze(listDatax[i_dat])
-        if vminList == None:
-            vmin = np.nanmin(data[:lim, :])
-        else:
-            vmin = vminList[i_dat]
-        if vmaxList == None:
-            vmax = np.nanmax(data[:lim, :])
-        else:
-            vmax = vmaxList[i_dat]
-        im = axs.imshow(
-            data[:lim, :],
-            cmap=cm.viridis,
-            interpolation="bicubic",
-            vmin=vmin,
-            vmax=vmax,
-            extent=[0, data_x[-1], tmax, 0],
-            aspect="auto",
-        )
-        divider = make_axes_locatable(axs)
-        cax = divider.append_axes("right", size="10%", pad=0.2)
-        cbar = fig.colorbar(im, cax=cax)
-        cbar.set_label(listCBLabel[i_dat])
-        ax = cbar.ax
-        text = ax.yaxis.label
-        font = matplotlib.font_manager.FontProperties(
-            family="times new roman", weight="bold", size=14
-        )
-        text.set_font_properties(font)
-        if listXAxisName is None:
-            if i_dat == 0:
-                pretty_labels("x", "t [s]", 12, listTitle[i_dat], ax=axs)
-            else:
-                pretty_labels("x", "", 12, listTitle[i_dat], ax=axs)
-        else:
-            if i_dat == 0:
-                pretty_labels(
-                    listXAxisName[i_dat], "t [s]", 12, listTitle[i_dat], ax=axs
-                )
-            else:
-                pretty_labels(
-                    listXAxisName[i_dat], "", 12, listTitle[i_dat], ax=axs
-                )
-
-        ax.set_xticks([])  # values
-        ax.set_xticklabels([])  # labels
-        if not i_dat == 0:
-            ax.set_yticks([])  # values
-            ax.set_yticklabels([])  # labels
-        for l in cbar.ax.yaxis.get_ticklabels():
-            l.set_weight("bold")
-            l.set_family("serif")
-            l.set_fontsize(12)
-    else:
-        for i_dat in range(len(listData)):
-            data = listData[i_dat]
-            data_x = np.squeeze(listDatax[i_dat])
-            if vminList == None:
-                vmin = np.nanmin(data[:lim, :])
-            else:
-                vmin = vminList[i_dat]
-            if vmaxList == None:
-                vmax = np.nanmax(data[:lim, :])
-            else:
-                vmax = vmaxList[i_dat]
-            im = axs[i_dat].imshow(
-                data[:lim, :],
-                cmap=cm.viridis,
-                interpolation="bicubic",
-                vmin=vmin,
-                vmax=vmax,
-                extent=[0, data_x[-1], tmax, 0],
-                aspect="auto",
-            )
-            divider = make_axes_locatable(axs[i_dat])
-            cax = divider.append_axes("right", size="10%", pad=0.2)
-            cbar = fig.colorbar(im, cax=cax)
-            cbar.set_label(listCBLabel[i_dat])
-            ax = cbar.ax
-            text = ax.yaxis.label
-            font = matplotlib.font_manager.FontProperties(
-                family="times new roman", weight="bold", size=14
-            )
-            text.set_font_properties(font)
-            if listXAxisName is None:
-                if i_dat == 0:
-                    pretty_labels(
-                        "x", "t [s]", 12, listTitle[i_dat], ax=axs[i_dat]
-                    )
-                else:
-                    pretty_labels("x", "", 12, listTitle[i_dat], ax=axs[i_dat])
-            else:
-                if i_dat == 0:
-                    pretty_labels(
-                        listXAxisName[i_dat],
-                        "t [s]",
-                        12,
-                        listTitle[i_dat],
-                        ax=axs[i_dat],
-                    )
-                else:
-                    pretty_labels(
-                        listXAxisName[i_dat],
-                        "",
-                        12,
-                        listTitle[i_dat],
-                        ax=axs[i_dat],
-                    )
-            axs[i_dat].set_xticks([])  # values
-            axs[i_dat].set_xticklabels([])  # labels
-            if not i_dat == 0:
-                axs[i_dat].set_yticks([])  # values
-                axs[i_dat].set_yticklabels([])  # labels
-            for l in cbar.ax.yaxis.get_ticklabels():
-                l.set_weight("bold")
-                l.set_family("serif")
-                l.set_fontsize(12)
-
-    if not globalTitle is None:
-        fig.tight_layout(rect=[0, 0.03, 1, 0.95])
-        plt.suptitle(
-            globalTitle,
-            fontsize=14,
-            fontweight="bold",
-            fontname="Times New Roman",
-        )
-
-
 def plotCollWeights(
     listDatax,
     listDatat,
diff --git a/requirements.txt b/requirements.txt
index c05a6a5..a08518e 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 tensorflow
-prettyPlot>=0.0.2
+prettyPlot>=0.0.5
 pandas<2.0
 corner
 tf2jax