diff --git a/Utilities/Python/fdsplotlib.py b/Utilities/Python/fdsplotlib.py
index ccd8cb695c1..bcaf62d6533 100644
--- a/Utilities/Python/fdsplotlib.py
+++ b/Utilities/Python/fdsplotlib.py
@@ -9,6 +9,7 @@
 Collection of functions for plotting and analysis
 """
 
+import os
 import sys
 import math
 import matplotlib
@@ -998,3 +999,145 @@ def dataplot(config_filename,**kwargs):
             print("Error in row {whichrow}, skipping case...".format(whichrow=irow+1))
             continue
 
+
+def add_version_string(ax, filename, plot_type='linear', scale_x=0.60, scale_y=1.02,
+                       font_name='Times', font_size=10):
+    """
+    Adds a version string to a matplotlib plot.
+
+    Parameters:
+    ax (matplotlib.axes.Axes): The axes to add the version string to.
+    filename (str): Path to the version string file.
+    plot_type (str): Type of plot ('loglog', 'semilogx', 'semilogy', or 'linear').
+    scale_x (float): Scaling factor for X-position.
+    scale_y (float): Scaling factor for Y-position.
+    font_name (str): Font name for the text.
+    font_interpreter (str): Interpreter type (not applicable in matplotlib, kept for compatibility).
+    font_size (int): Font size for the text.
+
+    """
+    if os.path.exists(filename):
+        with open(filename, 'r') as file:
+            version_str = file.read().strip()
+
+        x_lim = ax.get_xlim()
+        y_lim = ax.get_ylim()
+
+        if plot_type == 'loglog':
+            x_pos = 10**(np.log10(x_lim[0]) + scale_x * (np.log10(x_lim[1]) - np.log10(x_lim[0])))
+            y_pos = 10**(np.log10(y_lim[0]) + scale_y * (np.log10(y_lim[1]) - np.log10(y_lim[0])))
+        elif plot_type == 'semilogx':
+            x_pos = 10**(np.log10(x_lim[0]) + scale_x * (np.log10(x_lim[1]) - np.log10(x_lim[0])))
+            y_pos = y_lim[0] + scale_y * (y_lim[1] - y_lim[0])
+        elif plot_type == 'semilogy':
+            x_pos = x_lim[0] + scale_x * (x_lim[1] - x_lim[0])
+            y_pos = 10**(np.log10(y_lim[0]) + scale_y * (np.log10(y_lim[1]) - np.log10(y_lim[0])))
+        else:
+            x_pos = x_lim[0] + scale_x * (x_lim[1] - x_lim[0])
+            y_pos = y_lim[0] + scale_y * (y_lim[1] - y_lim[0])
+
+        ax.text(x_pos, y_pos, version_str, fontsize=font_size, fontname=font_name, verticalalignment='bottom')
+
+
+def get_plot_style(style="fds"):
+    """
+    Returns a dictionary of plot style parameters based on the specified style.
+
+    Parameters:
+    - style (str): The style to use ('fds', 'paper', etc.). Default is 'fds'.
+
+    Returns:
+    - dict: A dictionary containing plot style parameters.
+    """
+    if style == "fds":
+        return {
+            # Font properties
+            "Font_Name": "Times",
+            "Font_Interpreter": "TeX",
+            "Key_Font_Size": 12,
+            "Title_Font_Size": 16,
+            "Label_Font_Size": 16,
+            "Scat_Title_Font_Size": 14,
+            "Scat_Label_Font_Size": 14,
+            "Marker_Size": 4,
+            "D1_Marker_Size": 4,
+            "D2_Marker_Size": 4,
+
+            # Line properties
+            "Line_Width": 1.0,
+            "D1_Line_Width": 1.0,
+            "D2_Line_Width": 1.0,
+
+            # Plot properties
+            "Plot_Units": "inches",
+            "Plot_Width": 5.0,
+            "Plot_Height": 3.4,
+            "Plot_X": 1.2,
+            "Plot_Y": 0.8,
+            "Scat_Plot_Width": 4.75,
+            "Scat_Plot_Height": 4.75,
+            "Scat_Plot_X": 1.00,
+            "Scat_Plot_Y": 0.75,
+            "Subtitle_Text_Offset": 0.05,
+            "VerStr_Scale_X": 0.60,
+            "VerStr_Scale_Y": 1.05,
+
+            # Paper properties
+            "Paper_Units": "inches",
+            "Paper_Width": 6.5,
+            "Paper_Height": 4.5,
+            "Scat_Paper_Height": 6.0,
+            "Scat_Paper_Width": 6.0,
+
+            # Print properties
+            "Figure_Visibility": "off",
+            "Image_File_Type": "-dpdf",
+        }
+
+    elif style == "paper":
+        return {
+            # Font properties
+            "Font_Name": "Helvetica",
+            "Font_Interpreter": "TeX",
+            "Key_Font_Size": 16,
+            "Title_Font_Size": 20,
+            "Label_Font_Size": 20,
+            "Scat_Title_Font_Size": 14,
+            "Scat_Label_Font_Size": 14,
+            "Marker_Size": 10,
+            "D1_Marker_Size": 10,
+            "D2_Marker_Size": 10,
+
+            # Line properties
+            "Line_Width": 1.0,
+            "D1_Line_Width": 1.0,
+            "D2_Line_Width": 1.0,
+
+            # Plot properties
+            "Plot_Units": "normalized",
+            "Plot_X": 0.1500,
+            "Plot_Y": 0.1500,
+            "Plot_Width": 0.7750,
+            "Plot_Height": 0.8150 * 0.95,  # Adjusted for exponential y-axis tick labels
+            "Scat_Plot_Width": 4.75,
+            "Scat_Plot_Height": 4.75,
+            "Scat_Plot_X": 0.75,
+            "Scat_Plot_Y": 0.75,
+            "Subtitle_Text_Offset": 0.05,
+            "VerStr_Scale_X": 0.60,
+            "VerStr_Scale_Y": 1.05,
+
+            # Paper properties
+            "Paper_Units": "inches",
+            "Paper_Width": 8.0,
+            "Paper_Height": 6.0,
+            "Scat_Paper_Height": 6.0,
+            "Scat_Paper_Width": 6.0,
+
+            # Print properties
+            "Figure_Visibility": "on",
+            "Image_File_Type": "-dpdf",
+        }
+
+    else:
+        raise ValueError(f"Unknown style '{style}'. Please choose 'fds' or 'paper'.")
diff --git a/Utilities/Python/scripts/m2py_README.md b/Utilities/Python/scripts/m2py_README.md
new file mode 100644
index 00000000000..48d9f5c5a2b
--- /dev/null
+++ b/Utilities/Python/scripts/m2py_README.md
@@ -0,0 +1,38 @@
+# Matlab to Python conversion
+
+## 1. Use an AI code converter
+
+This step will get you 90% of the way there.  Then there usually a few extra steps to get a working script.
+
+## 2. Import `fdsplotlib`
+
+You will use this library for the Matlab equivalent of `plot_style` and `addverstr`.  At the top of your script add the following:
+
+```
+import fdsplotlib
+```
+
+For now, we will assume that you are running your scripts from the `fds/Utilities/Python/` directory where `fdsplotlib` lives.  If you need to run your python scripts from another directory, you need to add `fds/Utilities/Python/` to the path or else copy `fdsplotlibl.py` to your working directory.
+
+## 3. Add the plot style parameters
+
+Somewhere above where you make the plots add,
+
+```
+plot_style = fdsplotlib.get_plot_style("fds")
+plt.rcParams["font.family"] = plot_style["Font_Name"]
+plt.rcParams["font.size"] = plot_style["Label_Font_Size"]
+# print(plot_style) # uncomment to see list of parameters
+```
+
+This will polulate a Python "dictionary" called `plot_style` with all the parameters.  To see a list of the parameters uncomment the print statement.
+
+## 4. Add Git version string
+
+Similar to the Matlab plots, we need to add
+
+```
+fdsplotlib.add_version_string(ax, git_file, plot_type='linear')
+```
+
+where `ax` is the handle of the plot.
diff --git a/Utilities/Python/scripts/pyrolysis.py b/Utilities/Python/scripts/pyrolysis.py
index 99f9f71f680..b700f56e4c7 100644
--- a/Utilities/Python/scripts/pyrolysis.py
+++ b/Utilities/Python/scripts/pyrolysis.py
@@ -1,157 +1,141 @@
 import numpy as np
-import matplotlib.pyplot as plt
 from scipy.integrate import solve_ivp
+import matplotlib.pyplot as plt
 import os
-import sys
-
-# Add path to Verification/Pyrolysis
-sys.path.append('../../Verification/Pyrolysis')
-
-# Define global variables
-dTdt = None
-R0 = None
-E = None
-A = None
-residue = None
-
-def plot_style():
-    # Define plot style settings (placeholder)
-    plt.style.use('seaborn-darkgrid')
-    global Plot_Units, Plot_X, Plot_Y, Plot_Width, Plot_Height
-    global Font_Name, Label_Font_Size, Font_Interpreter
-    global Figure_Visibility, Paper_Units, Paper_Width, Paper_Height
-    Plot_Units = 'normalized'
-    Plot_X = 0.1
-    Plot_Y = 0.1
-    Plot_Width = 0.8
-    Plot_Height = 0.8
-    Font_Name = 'Arial'
-    Label_Font_Size = 12
-    Font_Interpreter = 'latex'
-    Figure_Visibility = 'on'
-    Paper_Units = 'inches'
-    Paper_Width = 6
-    Paper_Height = 4
-
-def addverstr(ax, filename, mode):
-    # Placeholder for addverstr function
-    # In Python, you might read the version string from the file and add text to the plot
-    if os.path.exists(filename):
-        with open(filename, 'r') as file:
-            ver_str = file.read().strip()
-            ax.text(0.05, 0.95, ver_str, transform=ax.transAxes, fontsize=8,
-                    verticalalignment='top')
-
-def reaction_rate(T, Y):
-    """
-    Defines the ODE system dY/dt = -A * Y * exp(-E / (R * T))
-    """
-    global A, E, R0
-    return -A * Y * np.exp(-E / (R0 * T))
-
-def main():
-    global dTdt, R0, E, A, residue
-
-    plt.close('all')
-    plt.figure()
-
-    plot_style()
+import pandas as pd
 
-    show_fds = True
+# include FDS plot styles, etc.
+import fdsplotlib
+import importlib
+importlib.reload(fdsplotlib) # use for development (while making changes to fdsplotlib.py)
 
-    for i_plot in range(1, 3):
-        fig, ax1 = plt.subplots()
-        fig.set_size_inches(Paper_Width, Paper_Height)
-        ax1.set_position([Plot_X, Plot_Y, Plot_Width, Plot_Height])
 
-        dTdt = 5.0 / 60.0  # degrees C per second
-        R0 = 8314.3  # J/(kmol*K)
-
-        if i_plot == 1:
-            n_components = 3
-            T_p = np.array([100.0 + 273.0, 300.0 + 273.0])
-            Y_0 = np.array([0.00, 1.00, 0.00])
-            delta_T = np.array([10.0, -80.0])
-            r_p = np.zeros(n_components - 1)
-            residue = np.array([0.0, 0.0])
+def reaction_rate(T, Y, dTdt, R0, E, A, residue):
+    """
+    Reaction rate function with safety checks and normalization of Y.
+    """
+    if T <= 0:
+        raise ValueError("Temperature T must be positive.")
+
+    # Ensure mass fractions are non-negative
+    Y = np.maximum(Y, 0.0)
+
+    r = np.zeros(len(Y))
+    r[0] = -A[0] * Y[0] * np.exp(-E[0] / (R0 * T)) / dTdt
+    r[1] = -A[1] * Y[1] * np.exp(-E[1] / (R0 * T)) / dTdt
+    r[2] = -residue[1] * r[1]
+
+    # Normalize Y to ensure the total remains valid
+    Y_sum = np.sum(Y)
+    if Y_sum > 0:
+        Y /= Y_sum
+
+    return r
+
+# Get plot style parameters
+plot_style = fdsplotlib.get_plot_style("fds")
+plt.rcParams["font.family"] = plot_style["Font_Name"]
+plt.rcParams["font.size"] = plot_style["Label_Font_Size"]
+
+# Define the base path for files
+base_path = "../../Verification/Pyrolysis/"  # Update this to the relative path of your files
+
+# Define parameters
+dTdt = 5.0 / 60.0
+R0 = 8314.3
+show_fds = True
+
+# Set up plots
+for i_plot in range(2):
+    fig, ax1 = plt.subplots(figsize=(plot_style["Paper_Width"], plot_style["Paper_Height"]))
+    ax2 = ax1.twinx()
+
+    # Define parameters for each plot
+    if i_plot == 0:
+        n_components = 3
+        T_p = np.array([100.0 + 273.0, 300.0 + 273.0])
+        Y_0 = np.array([0.00, 1.00, 0.00])
+        delta_T = np.array([10.0, -80.0])
+        r_p = np.array([0.0, 0.002])
+        residue = np.array([0.0, 0.0])
+    else:
+        n_components = 3
+        T_p = np.array([100.0 + 273.0, 300.0 + 273.0])
+        Y_0 = np.array([0.1, 0.9, 0.0])
+        delta_T = np.array([10.0, 80.0])
+        r_p = np.zeros(2)
+        residue = np.array([0.0, 0.2])
+
+    # Calculate A and E
+    A = np.zeros(n_components - 1)
+    E = np.zeros(n_components - 1)
+    for i in range(n_components - 1):
+        if delta_T[i] > 0:
+            r_p[i] = 2 * dTdt * (1.0 - residue[i]) / delta_T[i]
+        E[i] = np.exp(1.0) * r_p[i] * R0 * T_p[i] ** 2 / dTdt
+        A[i] = np.exp(1.0) * r_p[i] * np.exp(E[i] / (R0 * T_p[i]))
+
+    # Solve ODE
+    def odes(T, Y):
+        return reaction_rate(T, Y, dTdt, R0, E, A, residue)
+
+    sol = solve_ivp(odes, [273, 673], Y_0, method='BDF', rtol=1e-12, atol=1e-10)
+
+    T = sol.t
+    Y = sol.y.T
+    m = np.sum(Y, axis=1)
+
+    # Compute dm/dT
+    dmdT = np.gradient(m, T)
+
+    # Plot the solution
+    color1="tab:blue"
+    color2="tab:orange"
+    ax1.plot(T - 273, m, label="Normalized Mass", color=color1)
+    ax2.plot(T - 273, -dmdT * dTdt * 1000, label="Mass Loss Rate", color=color2)
+
+    # Load FDS solution
+    if show_fds:
+        if i_plot == 0:
+            fds_file = os.path.join(base_path, 'pyrolysis_1_devc.csv')
+            if not os.path.exists(fds_file):
+                print(f"Error: File {fds_file} does not exist. Skipping case.")
+                continue
         else:
-            n_components = 3
-            T_p = np.array([100.0 + 273.0, 300.0 + 273.0])
-            Y_0 = np.array([0.1, 0.9, 0.0])
-            delta_T = np.array([10.0, 80.0])
-            residue = np.array([0.0, 0.2])
-            r_p = np.zeros(n_components - 1)
-
-        # Calculate A and E
-        E = np.zeros(n_components - 1)
-        A = np.zeros(n_components - 1)
-        for i in range(n_components - 1):
-            if delta_T[i] > 0:
-                r_p[i] = 2.0 * dTdt * (1.0 - residue[i]) / delta_T[i]
-            E[i] = np.exp(1.0) * r_p[i] * R0 * T_p[i] ** 2 / dTdt
-            A[i] = np.exp(1.0) * r_p[i] * np.exp(E[i] / (R0 * T_p[i]))
-
-        # Solve the ODE dY/dT = f(T, Y)
-        # Since dY/dt = f(T, Y) and dT/dt is constant, dY/dT = f(T, Y) / dTdt
-        def dydT(T, Y):
-            return reaction_rate(T, Y) / dTdt
-
-        sol = solve_ivp(dydT, [273, 673], Y_0, method='BDF',
-                       rtol=1e-10, atol=1e-8)
-        T = sol.t
-        Y = sol.y.T
-
-        # Compute dm/dT
-        m = np.sum(Y, axis=1)
-        dmdT = np.gradient(m, T)
-        dmdT[0] = 0.0
-        dmdT[-1] = 0.0
-
-        # Plot the solution
-        TC = T - 273.0
-        ax2 = ax1.twinx()
-        line1, = ax1.plot(TC, m, label='Normalized Mass')
-        line2, = ax2.plot(TC, -dmdT * dTdt * 1000.0, label='Normalized Mass Loss Rate × 10³ (s⁻¹)', color='orange')
-        ax1.set_xlabel('Temperature (°C)', fontsize=Label_Font_Size, fontname=Font_Name, 
-                       fontfamily=Font_Interpreter)
-        ax1.set_ylabel('Normalized Mass', fontsize=Label_Font_Size, fontname=Font_Name, 
-                       fontfamily=Font_Interpreter)
-        ax2.set_ylabel('Normalized Mass Loss Rate × 10³ (s⁻¹)', fontsize=Label_Font_Size, fontname=Font_Name,
-                       fontfamily='LaTeX')
-        ax1.set_ylim([0, 1.1])
-        ax2.set_ylim([0, 2.2])
-        ax1.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
-        ax2.set_yticks([0, 0.4, 0.8, 1.2, 1.6, 2.0])
-        line1.set_linestyle('-')
-        line2.set_linestyle('-')
-
-        # Add the FDS solution
-        if show_fds:
-            filename = f'pyrolysis_{i_plot}_devc.csv'
-            if not os.path.exists(filename):
-                print(f'Error: File {filename} does not exist. Skipping case.')
+            fds_file = os.path.join(base_path, 'pyrolysis_2_devc.csv')
+            if not os.path.exists(fds_file):
+                print(f"Error: File {fds_file} does not exist. Skipping case.")
                 continue
-            FDS = np.genfromtxt(filename, delimiter=',', skip_header=2)
-            ax1.plot(FDS[:, 3], FDS[:, 1], 'b--', label='FDS Mass')
-            ax2.plot(FDS[:, 3], FDS[:, 2] * 500.0, 'r--', label='FDS Mass Loss Rate × 500')
 
-        # Add vertical lines to indicate the temperature peaks
-        for i in range(n_components - 1):
-            ax2.axvline(x=T_p[i] - 273.0, color='black', linewidth=1)
+        # Read the FDS data using pandas
+        fds_data = pd.read_csv(fds_file, skiprows=2, header=None)
+
+        # Convert to a NumPy array if required for further processing
+        fds_data = fds_data.to_numpy()
+
+    # Plot attributes
+    ax1.set_xlabel("Temperature (°C)",fontdict={"fontname": plot_style["Font_Name"], "fontsize": plot_style["Label_Font_Size"]})
+    ax1.set_ylabel("Normalized Mass",color=color1,fontdict={"fontname": plot_style["Font_Name"], "fontsize": plot_style["Label_Font_Size"]})
+    ax2.set_ylabel("Mass Loss Rate × 1000 (1/s)",color=color2,fontdict={"fontname": plot_style["Font_Name"], "fontsize": plot_style["Label_Font_Size"]})
+    ax1.tick_params(axis="y", colors=color1)
+    ax2.tick_params(axis="y", colors=color2)
+    ax1.set_ylim([0, 1.1])
+    ax2.set_ylim([0, 2.2])
 
-        # Add version string if file is available
-        chid = f'pyrolysis_{i_plot}'
-        git_filename = f'{chid}_git.txt'
-        addverstr(ax1, git_filename, 'linear')
+    # ax1.legend(loc="upper left")
+    # ax2.legend(loc="upper right")
 
-        # Create the PDF files
-        fig.set_visible(True if Figure_Visibility == 'on' else False)
-        fig.set_size_inches(Paper_Width * 1.1, Paper_Height)
-        pdf_filename = f'../../Manuals/FDS_User_Guide/SCRIPT_FIGURES/{chid}.pdf'
-        plt.savefig(pdf_filename, format='pdf')
+    # Add vertical lines for temperature peaks
+    for i in range(n_components - 1):
+        ax2.axvline(T_p[i] - 273, color="black", linestyle="-", linewidth=1)
 
-    plt.show()
+    # Add version sting
+    chid = 'pyrolysis_1' if i_plot == 0 else 'pyrolysis_2'
+    git_file = os.path.join(base_path, f"{chid}_git.txt")
+    fdsplotlib.add_version_string(ax1, git_file, plot_type='linear')
 
-if __name__ == "__main__":
-    main()
+    plt.tight_layout()
 
+    # Save the figure
+    plt.savefig(f"../../Manuals/FDS_User_Guide/SCRIPT_FIGURES/pyrolysis_{i_plot + 1}.pdf")
+    plt.close()
diff --git a/Utilities/Python/scripts/pyrolysis_notes.txt b/Utilities/Python/scripts/pyrolysis_notes.txt
deleted file mode 100644
index 23c999121a1..00000000000
--- a/Utilities/Python/scripts/pyrolysis_notes.txt
+++ /dev/null
@@ -1,20 +0,0 @@
-
-Functions plot_style and addverstr: In the MATLAB code, these are likely custom functions for setting plot styles and adding version strings. In the Python translation, plot_style initializes plot parameters using matplotlib, and addverstr is a placeholder that reads a version string from a file and adds it to the plot. You may need to implement the actual logic based on your specific requirements.
-
-Global Variables: MATLAB's global variables are translated to module-level variables in Python. However, it's generally recommended to pass variables explicitly to functions to avoid potential issues with state management.
-
-ODE Solver: MATLAB's ode15s is translated to scipy.integrate.solve_ivp with the BDF method, which is suitable for stiff ODEs. The ODE function reaction_rate is adjusted to account for the difference in defining dY/dT based on the temperature rate dTdt.
-
-File Reading: MATLAB's csvread is replaced with numpy.genfromtxt, which is more flexible and robust for reading CSV files in Python.
-
-Plotting: MATLAB's plotyy is handled using Matplotlib's twinx to create a secondary y-axis. Labels, legends, and styles are set using Matplotlib's API.
-
-Indexing: MATLAB indices start at 1, whereas Python indices start at 0. Adjustments have been made to account for this difference, especially in loop ranges and array indexing.
-
-Error Handling: The MATLAB display and return statements are translated to Python's print and continue statements to handle missing files gracefully without terminating the entire script.
-
-Units and Calculations: Ensure that the units (e.g., temperature in Kelvin vs. Celsius) are consistently handled throughout the translation to maintain the integrity of the computations.
-
-Placeholder Implementations: Certain functions like plot_style and addverstr are provided as placeholders. You'll need to implement their specific functionalities based on your original MATLAB definitions or requirements.
-
-File Paths: The addpath in MATLAB is translated to modifying sys.path in Python. Ensure that the relative paths are correct based on your project's directory structure.