diff --git a/Project.toml b/Project.toml
index ebaae8a..ec73e09 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ElectrochemicalKinetics"
 uuid = "a2c6e634-85ca-418a-9c67-9b5417ce2d04"
 authors = ["Rachel Kurchin <rkurchin@cmu.edu>", "Holden Parks <hparks@andrew.cmu.edu>", "Dhairya Gandhi <dhairya@juliacomputing.com>"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
diff --git a/src/phase_diagrams.jl b/src/phase_diagrams.jl
index 5acbe60..c082408 100644
--- a/src/phase_diagrams.jl
+++ b/src/phase_diagrams.jl
@@ -63,13 +63,21 @@ function find_phase_boundaries(I, km::KineticModel; guess=[0.05, 0.95], intercal
     x1.zero
 end
 
-# TODO: add T_max to handle low-T BV cases
-function phase_diagram(km::KineticModel; I_step=10, verbose=false, intercalate=true, kwargs...)
+"""
+    phase_diagram(km; I_step=1, I_max=Inf, verbose=false, intercalate=true, kwargs...)
+
+Construct electrochemical phase diagram for the model `km` as a function of composition and current, in steps of `I_step`.
+    
+NOTE 1: appropriate values of `I_step` depend strongly on the prefactor of your model. For example, for `ButlerVolmer` with a prefactor of 1, the phase diagram at T=330K closes at I=5, but with a prefactor of 10, it reaches up to I=44. 
+
+NOTE 2: at lower temperatures (<=320K or so), ButlerVolmer models with the default thermodynamic parameters have a two-phase region at every current, so setting a finite value of I_max is necessary for this function to finish running.
+"""
+function phase_diagram(km::KineticModel; I_step=1, I_max=Inf, verbose=false, intercalate=true, kwargs...)
     I = 0
     pbs_here = find_phase_boundaries(I, km; intercalate=intercalate, kwargs...)
     pbs = pbs_here'
     I_vals = [0]
-    while abs(pbs_here[2] - pbs_here[1]) > 1e-3
+    while abs(pbs_here[2] - pbs_here[1]) > 1e-3 && I < I_max
         I = I + I_step
         if verbose
             println("Solving at I=", I, "...")