diff --git a/examples/1D_convergence/README.md b/examples/1D_convergence/README.md
new file mode 100644
index 0000000000..3b9963e048
--- /dev/null
+++ b/examples/1D_convergence/README.md
@@ -0,0 +1,11 @@
+This example case contains an automated convergence test using a 1D, two-component advection case.
+The case can be run by executing the bash script `./submitJobs.sh` in a terminal after enabling
+execution permissions with `chmod +x ./submitJobs.sh` and setting the `ROOT_DIR` and `MFC_DIR`
+variables. By default the script runs the case for 6 different grid resolutions with 1st,
+3rd, and 5th, order spatial reconstructions. These settings can be modified by editing the variables
+at the top of the script. You can also run different model equations by setting the `ME` variable
+and different Riemann solvers by setting the `RS` variable.
+
+Once the simulations have been run, you can generate convergence plots with matplotlib by running
+`python3 plot.py` in a terminal. This will generate plots of the L1, L2, and Linf error norms and
+save the results to `errors.csv`.
diff --git a/examples/1D_convergence/case.py b/examples/1D_convergence/case.py
new file mode 100755
index 0000000000..fe1ed3c976
--- /dev/null
+++ b/examples/1D_convergence/case.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+import math
+import json
+import argparse
+
+# Parsing command line arguments
+parser = argparse.ArgumentParser(description="Generate JSON case configuration for two-fluid convergence simulation.")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC's toolchain's internal state.")
+parser.add_argument("--order", type=int, default=5, help="WENO order (default: 5)")
+parser.add_argument("--meqns", type=int, default=2, help="Model equations (default: 2 (five-equation model))")
+parser.add_argument("--rs", type=int, default=2, help="Riemann solver (default: 2 (HLLC))")
+parser.add_argument("-N", type=int, default=1024, help="Number of grid points (default: 1024)")
+
+args = parser.parse_args()
+
+# Numerical setup
+Nx = args.N - 1
+dx = 1.0 / (1.0 * (Nx + 1))
+
+Tend, Nt = 5.0, 200000
+mydt = Tend / (1.0 * Nt)
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0e00,
+            "x_domain%end": 1.0e00,
+            "m": Nx,
+            "n": 0,
+            "p": 0,
+            "dt": mydt,
+            "t_step_start": 0,
+            "t_step_stop": int(Nt),
+            "t_step_save": int(math.ceil(Nt / 10.0)),
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": args.meqns,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "weno_order": args.order,
+            "weno_eps": dx**2,
+            "weno_Re_flux": "F",
+            "weno_avg": "F",
+            "mapped_weno": "F",
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": args.rs,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1 L
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%vel(1)": 1.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": f"0.5 - 0.5*sin(2*pi*x)",
+            "patch_icpp(1)%alpha(1)": f"0.5 - 0.5*sin(2*pi*x)",
+            "patch_icpp(1)%alpha_rho(2)": f"0.5 + 0.5*sin(2*pi*x)",
+            "patch_icpp(1)%alpha(2)": f"0.5 + 0.5*sin(2*pi*x)",
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(2)%gamma": 1.0e00 / (1.4 - 1.0e00),
+            "fluid_pp(2)%pi_inf": 0.0,
+        }
+    )
+)
diff --git a/examples/1D_convergence/plot.py b/examples/1D_convergence/plot.py
new file mode 100644
index 0000000000..60747b6965
--- /dev/null
+++ b/examples/1D_convergence/plot.py
@@ -0,0 +1,69 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+N = [32, 64, 128, 256, 512, 1024]
+Ord = [1, 3, 5]
+
+errors = np.nan * np.zeros((len(N), len(Ord), 3))
+
+TEND = 200000
+
+for i in range(len(N)):
+    for j in range(len(Ord)):
+
+        sim_a1 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.5.00.{TEND}.dat", sep=r"\s+", header=None, names=["x", "y"])
+        sim_a2 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.6.00.{TEND}.dat", sep=r"\s+", header=None, names=["x", "y"])
+
+        exact_a1 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.5.00.000000.dat", sep=r"\s+", header=None, names=["x", "y"])
+        exact_a2 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.6.00.000000.dat", sep=r"\s+", header=None, names=["x", "y"])
+
+        ## 2 norm
+        errors[i, j, 0] = np.linalg.norm(sim_a1.y - exact_a1.y) / np.sqrt(N[i])
+        errors[i, j, 0] += np.linalg.norm(sim_a2.y - exact_a2.y) / np.sqrt(N[i])
+
+        ## 1 norm
+        errors[i, j, 1] = 1 / N[i] * np.sum(np.abs(sim_a1.y - exact_a1.y))
+        errors[i, j, 1] += 1 / N[i] * np.sum(np.abs(sim_a2.y - exact_a2.y))
+
+        ## Inf norm
+        errors[i, j, 2] = np.max([np.nanmax(np.abs(sim_a1.y - exact_a1.y)), np.nanmax(np.abs(sim_a2.y - exact_a2.y))])
+
+fig, ax = plt.subplots(1, 3, figsize=(12, 8), sharex=True)
+
+colors = ["blue", "green", "red", "purple"]
+
+ref = np.nan * np.zeros((len(N), len(Ord)))
+
+for i in range(3):
+    ax[i].plot(N, 30 / np.array(N) ** 1, label="Slope = -1", color=colors[0])
+    ax[i].plot(N, 3000 / np.array(N) ** 3, label="Slope = -3", color=colors[1])
+    ax[i].plot(N, 5000 / np.array(N) ** 5, label="Slope = -5", color=colors[2])
+
+for j in range(len(Ord)):
+    ax[0].plot(N, errors[:, j, 0], "o-", color=colors[j])
+    ax[0].set_xscale("log", base=2)
+    ax[0].set_yscale("log")
+    ax[0].set_title("||error||_2")
+    ax[0].legend()
+
+    ax[1].plot(N, errors[:, j, 1], "o-", color=colors[j])
+    ax[1].set_xscale("log", base=2)
+    ax[1].set_yscale("log")
+    ax[1].set_title("||error||_1")
+
+    ax[2].plot(N, errors[:, j, 2], "o-", color=colors[j])
+    ax[2].set_xscale("log", base=2)
+    ax[2].set_yscale("log")
+    ax[2].set_title("||error||_inf")
+
+plt.tight_layout()
+plt.show()
+
+errors = np.column_stack((N, errors[:, :, 0]))
+errors = np.column_stack((errors, 7 / np.array(N) ** 1))
+errors = np.column_stack((errors, 700 / np.array(N) ** 3))
+errors = np.column_stack((errors, 2000 / np.array(N) ** 5))
+
+df = pd.DataFrame(errors, columns=["N", "1", "3", "5", "R1", "R3", "R5"], index=N)
+df.to_csv("errors.csv", index=False)
diff --git a/examples/1D_convergence/submitJobs.sh b/examples/1D_convergence/submitJobs.sh
new file mode 100755
index 0000000000..0f5315fa9e
--- /dev/null
+++ b/examples/1D_convergence/submitJobs.sh
@@ -0,0 +1,28 @@
+#!/usr/bin/env sh
+Nx=(32 64 128 256 512 1024)
+Order=(1 3 5)
+
+ME=2 # Model equations = 2 for five-equation model
+RS=2 # Riemann solver = 2 for HLLC
+
+# Get the directory of the script itself
+ROOT_DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
+# Assume the script is in examples/1D_convergence, so MFC_DIR is two levels up
+MFC_DIR="$(dirname "$(dirname "$ROOT_DIR")")"
+
+for i in "${Nx[@]}"; do
+    for j in "${Order[@]}"; do
+        rm -rf N${i}_O${j}
+        mkdir N${i}_O${j}
+        cp case.py N${i}_O${j}/
+    done
+done
+
+cd $MFC_DIR
+
+for i in "${Nx[@]}"; do
+    for j in "${Order[@]}"; do
+        ./mfc.sh run $ROOT_DIR/N${i}_O${j}/case.py --case-optimization --no-debug -- --order $j -N $i --meqns $ME --rs $RS
+    done
+done
+
diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py
index fb30a21f61..269c8438bf 100644
--- a/toolchain/mfc/case.py
+++ b/toolchain/mfc/case.py
@@ -142,7 +142,7 @@ def rhs_replace(match):
                     'r':     f'patch_icpp({pid})%radius',  'eps':   f'patch_icpp({pid})%epsilon', 'beta':  f'patch_icpp({pid})%beta',
                     'tau_e': f'patch_icpp({pid})%tau_e',   'radii': f'patch_icpp({pid})%radii',
 
-                    'e' : f'{math.e}', 'pi': f'{math.pi}',
+                    'e' : f'{math.e}',
                 }.get(match.group(), match.group())
 
             lines = []
diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py
index 6fb00781be..fa23933b69 100644
--- a/toolchain/mfc/test/cases.py
+++ b/toolchain/mfc/test/cases.py
@@ -1021,7 +1021,8 @@ def foreach_example():
                            "3D_TaylorGreenVortex_analytical",
                            "3D_IGR_TaylorGreenVortex_nvidia",
                            "2D_backward_facing_step",
-                           "2D_forward_facing_step"]
+                           "2D_forward_facing_step",
+                           "1D_convergence"]
             if path in casesToSkip:
                 continue
             name = f"{path.split('_')[0]} -> Example -> {'_'.join(path.split('_')[1:])}"