diff --git a/Manifest.toml b/Manifest.toml
index 3a279f8..4e4a3fb 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "307cac97db0b7035e7716e25208fc63e328a2438"
+project_hash = "ce9748a73286acaa3182365f1a6814016433d441"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
diff --git a/Project.toml b/Project.toml
index e4fb426..cafde87 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 IterableTables = "1c8ee90f-4401-5389-894e-7a04a3dc0f4d"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
diff --git a/scripts/runcode.jl b/scripts/runcode.jl
index d049bfb..d4fdeb9 100644
--- a/scripts/runcode.jl
+++ b/scripts/runcode.jl
@@ -268,10 +268,12 @@ julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialCo
         df = DataFrame(sol)
         write_csv(df, erupt_saved, path)
         if plotfig
-            plot_figs(df, path)
+            title = latexstringtitle(composition, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth)
+            plot_combined_fig(df, path, title)
+            plot_dual_axis(df, param_saved_var.storeTime, param_saved_var.storeEps_x, path)
         end
-        write_ϵx_csv(param_saved_var.storeTime, param_saved_var.storeEps_x, path)
-        plot_ϵx(param_saved_var.storeTime, param_saved_var.storeEps_x, path)
+        # write_ϵx_csv(param_saved_var.storeTime, param_saved_var.storeEps_x, path)
+        # plot_ϵx(param_saved_var.storeTime, param_saved_var.storeEps_x, path)
         return df
     end
     return df
diff --git a/src/Chamber.jl b/src/Chamber.jl
index bbd7499..b0be27e 100644
--- a/src/Chamber.jl
+++ b/src/Chamber.jl
@@ -11,6 +11,7 @@ using Roots
 using SpecialFunctions
 using Sundials
 using TimerOutputs
+using LaTeXStrings
 using Logging
 include("utils.jl")
 include("./glq_points_weights.jl")
@@ -135,6 +136,10 @@ export @memoize,
     build_rho_rc,
     write_csv,
     plot_figs,
+    plot_sol_year_unit,
+    latexstringtitle,
+    plot_combined_fig,
+    plot_dual_axis,
     plot_ϵx,
     write_ϵx_csv,
     check_for_duplicates,
diff --git a/src/ic_finder.jl b/src/ic_finder.jl
index d8ac0aa..f7999c5 100644
--- a/src/ic_finder.jl
+++ b/src/ic_finder.jl
@@ -268,7 +268,7 @@ function IC_Finder(
         end
     end
 
-    if eps_g0 <= 0 || X_co20 < 0
+    if eps_g0 <= 0 || X_co20 < 0 || any(isnan, [X_co20, eps_g0])
         X_co20 = 0.0
         eps_g0 = 0.0
         mco2_diss = m_co2_melt
diff --git a/src/plot_figs.jl b/src/plot_figs.jl
index 39acc06..d4ded76 100644
--- a/src/plot_figs.jl
+++ b/src/plot_figs.jl
@@ -4,12 +4,38 @@ function plot_sol(df, x_axis_col, y_axis_col)
         df[:, y_axis_col];
         xaxis=x_axis_col,
         yaxis=y_axis_col,
-        label="Chamber.jl",
         linewidth=2,
         marker=(:x, 3, Plots.stroke(2)),
     )
 end
 
+function plot_sol_year_unit(df, xlabel, ylabel)
+    # Units of outputs: Time -> sec, Presure -> Pa, Volume -> m³
+    divisors = Dict(
+        "time" => 31536000,
+        "P+dP" => 1e6,
+        "V" => 1e9,
+        "total_mass" => 1e12,
+    )
+    labels = Dict(
+        "time" => "Time (yr)",
+        "P+dP" => "Pressure (MPa)",
+        "T" => "Temperature (K)",
+        "V" => "Volume (km³)",
+        "total_mass" => "Total mass (×10¹² kg)",
+    )
+    convert_unit(label) = df[:, label] / get(divisors, label, 1.0)
+    getaxis(label) = get(labels, label, label)
+
+    return plot(
+        convert_unit(xlabel), convert_unit(ylabel);
+        xaxis=getaxis(xlabel),
+        yaxis=getaxis(ylabel),
+        legend=false,
+        linewidth=2,
+    )
+end
+
 function plot_figs(df::DataFrame, path::String)::Nothing
     ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
     savefig′(label, path) = savefig(plot_sol(df, "time", label), path)
@@ -19,10 +45,19 @@ function plot_figs(df::DataFrame, path::String)::Nothing
     return nothing
 end
 
-function plot_ϵx(storeTime::Vector{Float64}, storeEps_g::Vector{Float64}, path::String)::Nothing
+function plot_combined_fig(df::DataFrame, path::String, title::LaTeXString)::Nothing
+    ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
+    plot_names = ["P+dP", "T", "eps_g", "V", "X_CO2", "total_mass"]
+    plts = [plot_sol_year_unit(df, "time", l) for l in plot_names]
+    combined_plot = plot(plts..., layout=(2, 3), size=(900, 500), dpi=200, margin = 3Plots.mm, plot_title=title, plot_titlefontsize=13)
+    savefig(combined_plot, joinpath(path, "combined_plots.png"))
+    return nothing
+end
+
+function plot_ϵx(storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing
     ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
     savefig(plot(
-            storeTime, storeEps_g;
+            storeTime, storeEps_x;
             xaxis="time",
             yaxis="eps_x",
             label="Chamber.jl",
@@ -31,3 +66,62 @@ function plot_ϵx(storeTime::Vector{Float64}, storeEps_g::Vector{Float64}, path:
         "$path/eps_x.png")
     return nothing
 end
+
+function plot_dual_axis(df::DataFrame, storeTime::Vector{Float64}, storeEps_x::Vector{Float64}, path::String)::Nothing
+    color1 = RGB(20/255, 125/255, 204/255)
+    color2 = RGB(255/255, 100/255, 0/255)
+    ENV["GKSwstype"] = "100"  # magic environmental variable for Plots
+
+    common_opt1 = (
+        xaxis="Time (yr)",
+        linewidth=3,
+        legend=:false,
+        guidefontsize=18,
+        tickfontsize=12,
+        linecolor=color1,
+        fg_color_guide=color1,
+        fg_color_text=color1,
+    )
+
+    common_opt2 = (
+        linewidth=3,
+        legend=:false,
+        guidefontsize=18,
+        tickfontsize=12,
+        linecolor=color2,
+        fg_color_guide=color2,
+        fg_color_text=color2,
+    )
+
+    plot(
+        df[!, 1]/31536000,
+        df[!, 2]/1e6;
+        common_opt1...,
+        yaxis="Pressure (MPa)",
+        )
+    p1 = plot!(
+        twinx(),
+        df[!, 1]/31536000,
+        df[!, 3];
+        common_opt2...,
+        yaxis="Temperature (K)",
+    )
+
+    plot(
+        df[!, 1]/31536000,
+        df[!, 4];
+        common_opt1...,
+        yaxis=L"𝜺_g",
+        )
+    p2 = plot!(
+        twinx(),
+        storeTime/31536000,
+        storeEps_x;
+        common_opt2...,
+        yaxis=L"𝜺_X",
+    )
+
+    savefig(p1, joinpath(path, "P_T.png"))
+    savefig(p2, joinpath(path, "eps_x_eps_g.png"))
+    return nothing
+end
diff --git a/src/utils.jl b/src/utils.jl
index d0f3a97..d82d1ad 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -506,3 +506,5 @@ function check_for_duplicates(
     )
     return nothing
 end
+
+latexstringtitle(composition, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth) = L"\textrm{%$(string(typeof(composition)))\ \ H_2O:%$(InitialConc_H2O*100)\ \%,\ CO_2: %$(InitialConc_CO2*1000000)\ ppm,\ depth:%$(depth/1000)\ km,\ vfr: 10^{%$(log_vfr)}\ km^{3}/yr,\ V_0:10^{%$(log_volume_km3)}\ km^{3}}"