diff --git a/src/simulation.py b/src/simulation.py
index 966a4a9..c46f9d3 100644
--- a/src/simulation.py
+++ b/src/simulation.py
@@ -149,11 +149,9 @@ def run_single_simulation(args):
     except ValueError:  # Loss rate too high
         return None
 
-def simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, lambda_vals):
+def simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, lambda_vals, ax_rt):
     C_values = [1, 2, 3, 6]
 
-    plt.figure(figsize=(12, 8))
-
     mean_at_lambda1 = {}
     loss_data = {}
 
@@ -205,8 +203,8 @@ def simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, l
                     print(f"Stopped at λ={lambda_val:.2f} - no successful run")
                     break
 
-            plt.plot(lambda_points, means, label=f'C={c}')
-            plt.fill_between(lambda_points, ci_lower, ci_upper, alpha=0.2)
+            ax_rt.plot(lambda_points, means, label=f'C={c}')
+            ax_rt.fill_between(lambda_points, ci_lower, ci_upper, alpha=0.2)
 
             # store response time for lamba = 1
             if 1 in lambda_points:
@@ -229,17 +227,18 @@ def simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, l
 
     # plot curves
     if len(lambda_vals)>1:
-        plt.xlabel('Arrival Rate (λ)')
-        plt.ylabel('Mean Response Time')
-        plt.title(f'Mean Response Time vs Arrival Rate ({num_runs} runs, 95% CI)')
-        plt.legend()
-        plt.grid(True)
+        ax_rt.set_xlim(left=0)
+        ax_rt.set_xlabel('Arrival Rate (λ)')
+        ax_rt.set_ylabel('Mean Response Time')
+        ax_rt.set_title(f'Mean Response Time vs Arrival Rate ({num_runs} runs, {int(confidence_level*100)}% CI)')
+        ax_rt.legend()
+        ax_rt.grid(True)
 
     return loss_data
 
 def simulation_wrapper2(simulation_time, num_runs, min_runs, confidence_level, lambda_vals, max_loss_value=0.1):
-    loss_data = simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, lambda_vals)
-    plt.show()
+    fig, (ax_rt, ax_loss) = plt.subplots(2,1, figsize=(12,10), gridspec_kw={'height_ratios': [4, 3], 'hspace': 0.4}, sharex=False)
+    loss_data = simulation_wrapper1(simulation_time, num_runs, min_runs, confidence_level, lambda_vals, ax_rt)
 
     for c, (lams, losses) in loss_data.items():
         lambda_vals2 = []
@@ -300,14 +299,16 @@ def simulation_wrapper2(simulation_time, num_runs, min_runs, confidence_level, l
                     print(f"Stopped at λ={lambda_val:.4f} - no successful run")
                     break
 
-            plt.plot(lambda_points, losses, label=f'C={c}')
-            plt.fill_between(lambda_points, ci_lower, ci_upper, alpha=0.2)
+            ax_loss.plot(lambda_points, losses, label=f'C={c}')
+            ax_loss.fill_between(lambda_points, ci_lower, ci_upper, alpha=0.2)
+
+    ax_loss.set_xlim(left=2)
+    ax_loss.set_xlabel('Arrival Rate (λ)')
+    ax_loss.set_ylabel('Mean Loss Rate')
+    ax_loss.set_title(f'Mean Loss Rate vs Arrival Rate ({num_runs} runs, {int(confidence_level*100)}% CI)')
+    ax_loss.legend()
+    ax_loss.grid(True)
 
-    plt.xlabel('Arrival Rate (λ)')
-    plt.ylabel('Mean Loss Rate')
-    plt.title(f'Mean Loss Rate vs Arrival Rate ({num_runs} runs, 95% CI)')
-    plt.legend()
-    plt.grid(True)
     plt.show()