diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
new file mode 100644
index 0000000..030c5a4
--- /dev/null
+++ b/.github/workflows/Documentation.yml
@@ -0,0 +1,18 @@
+name: "Documentation"
+
+on:
+  push:
+    branches:
+      - master
+    tags: '*'
+  pull_request:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ github.ref_name != github.event.repository.default_branch || github.ref != 'refs/tags/v*' }}
+
+jobs:
+  build-and-deploy-docs:
+    name: "Documentation"
+    uses: "SciML/.github/.github/workflows/documentation.yml@v1"
+    secrets: "inherit"
diff --git a/docs/Project.toml b/docs/Project.toml
index 6118f9b..57f966a 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,11 @@
 [deps]
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NeuralLyapunov = "61252e1c-87f0-426a-93a7-3cdb1ed5a867"
+NeuralPDE = "315f7962-48a3-4962-8226-d0f33b1235f0"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/make.jl b/docs/make.jl
index 5779690..5c82102 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -4,6 +4,22 @@ using Documenter
 DocMeta.setdocmeta!(
     NeuralLyapunov, :DocTestSetup, :(using NeuralLyapunov); recursive = true)
 
+MANUAL_PAGES = [
+    "man.md",
+    "man/pdesystem.md",
+    "man/minimization.md",
+    "man/decrease.md",
+    "man/structure.md",
+    "man/roa.md",
+    "man/policy_search.md",
+    "man/local_lyapunov.md"
+]
+DEMONSTRATION_PAGES = [
+    "demos/damped_SHO.md",
+    "demos/roa_estimation.md",
+    "demos/policy_search.md"
+]
+
 makedocs(;
     modules = [NeuralLyapunov],
     authors = "Nicholas Klugman <13633349+nicholaskl97@users.noreply.github.com> and contributors",
@@ -14,7 +30,9 @@ makedocs(;
         assets = String[]
     ),
     pages = [
-        "Home" => "index.md"
+        "Home" => "index.md",
+        "Manual" => vcat(MANUAL_PAGES, hide("man/internals.md")),
+        "Demonstrations" => DEMONSTRATION_PAGES
     ]
 )
 
diff --git a/docs/src/demos/damped_SHO.md b/docs/src/demos/damped_SHO.md
new file mode 100644
index 0000000..46f7104
--- /dev/null
+++ b/docs/src/demos/damped_SHO.md
@@ -0,0 +1,282 @@
+# Damped Simple Harmonic Oscillator
+
+Let's train a neural network to prove the exponential stability of the damped simple harmonic oscillator (SHO).
+
+The damped SHO is represented by the system of differential equations
+```math
+\begin{align*}
+    \frac{dx}{dt} &= v \\
+    \frac{dv}{dt} &= -2 \zeta \omega_0 v - \omega_0^2 x
+\end{align*}
+```
+where ``x`` is the position, ``v`` is the velocity, ``t`` is time, and ``\zeta, \omega_0`` are parameters.
+
+We'll consider just the box domain ``x \in [-5, 5], v \in [-2, 2]``.
+
+## Copy-Pastable Code
+
+```julia
+using NeuralPDE, Lux, NeuralLyapunov
+using Optimization, OptimizationOptimisers, OptimizationOptimJL
+using Random
+
+Random.seed!(200)
+
+######################### Define dynamics and domain ##########################
+
+"Simple Harmonic Oscillator Dynamics"
+function f(state, p, t)
+    ζ, ω_0 = p
+    pos = state[1]
+    vel = state[2]
+    vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos)
+end
+lb = [-5.0, -2.0];
+ub = [ 5.0,  2.0];
+p = [0.5, 1.0];
+fixed_point = [0.0, 0.0];
+dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0]))
+
+####################### Specify neural Lyapunov problem #######################
+
+# Define neural network discretization
+dim_state = length(lb)
+dim_hidden = 10
+dim_output = 3
+chain = [Lux.Chain(
+             Dense(dim_state, dim_hidden, tanh),
+             Dense(dim_hidden, dim_hidden, tanh),
+             Dense(dim_hidden, 1)
+         ) for _ in 1:dim_output]
+
+# Define training strategy
+strategy = QuasiRandomTraining(1000)
+discretization = PhysicsInformedNN(chain, strategy)
+
+# Define neural Lyapunov structure
+structure = NonnegativeNeuralLyapunov(
+    dim_output;
+    δ = 1e-6
+)
+minimization_condition = DontCheckNonnegativity(check_fixed_point = true)
+
+# Define Lyapunov decrease condition
+# Damped SHO has exponential decrease at a rate of k = ζ * ω_0, so we train to certify that
+decrease_condition = ExponentialDecrease(prod(p))
+
+# Construct neural Lyapunov specification
+spec = NeuralLyapunovSpecification(
+    structure,
+    minimization_condition,
+    decrease_condition
+)
+
+############################# Construct PDESystem #############################
+
+@named pde_system = NeuralLyapunovPDESystem(
+    dynamics,
+    lb,
+    ub,
+    spec;
+    p = p
+)
+
+######################## Construct OptimizationProblem ########################
+
+prob = discretize(pde_system, discretization)
+
+########################## Solve OptimizationProblem ##########################
+
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500)
+prob = Optimization.remake(prob, u0 = res.u)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 500)
+
+###################### Get numerical numerical functions ######################
+net = discretization.phi
+θ = res.u.depvar
+
+V, V̇ = get_numerical_lyapunov_function(
+    net,
+    θ,
+    structure,
+    f,
+    fixed_point;
+    p = p
+)
+```
+
+## Detailed description
+
+In this example, we set the dynamics up as an `ODEFunction` and use a `SciMLBase.SymbolCache` to tell the ultimate `PDESystem` what to call our state and parameter variables.
+
+```@setup SHO
+using Random
+
+Random.seed!(200)
+```
+
+```@example SHO
+using NeuralPDE # for ODEFunction and SciMLBase.SymbolCache
+
+"Simple Harmonic Oscillator Dynamics"
+function f(state, p, t)
+    ζ, ω_0 = p
+    pos = state[1]
+    vel = state[2]
+    vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos)
+end
+lb = [-5.0, -2.0];
+ub = [ 5.0,  2.0];
+p = [0.5, 1.0];
+fixed_point = [0.0, 0.0];
+dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0]))
+```
+
+Setting up the neural network using Lux and NeuralPDE training strategy is no different from any other physics-informed neural network problem.
+For more on that aspect, see the [NeuralPDE documentation](https://docs.sciml.ai/NeuralPDE/stable/).
+
+```@example SHO
+using Lux
+
+# Define neural network discretization
+dim_state = length(lb)
+dim_hidden = 10
+dim_output = 3
+chain = [Lux.Chain(
+             Dense(dim_state, dim_hidden, tanh),
+             Dense(dim_hidden, dim_hidden, tanh),
+             Dense(dim_hidden, 1)
+         ) for _ in 1:dim_output]
+```
+
+```@example SHO
+# Define training strategy
+strategy = QuasiRandomTraining(1000)
+discretization = PhysicsInformedNN(chain, strategy)
+```
+
+We now define our Lyapunov candidate structure along with the form of the Lyapunov conditions we'll be using.
+
+For this example, let's use a Lyapunov candidate
+```math
+V(x) = \lVert \phi(x) \rVert^2 + \delta \log \left( 1 + \lVert x \rVert^2 \right),
+```
+which structurally enforces nonnegativity, but doesn't guarantee ``V([0, 0]) = 0``.
+We therefore don't need a term in the loss function enforcing ``V(x) > 0 \, \forall x \ne 0``, but we do need something enforcing ``V([0, 0]) = 0``.
+So, we use [`DontCheckNonnegativity(check_fixed_point = true)`](@ref).
+
+To train for exponential decrease we use [`ExponentialDecrease`](@ref), but we must specify the rate of exponential decrease, which we know in this case to be ``\zeta \omega_0``.
+
+```@example SHO
+using NeuralLyapunov
+
+# Define neural Lyapunov structure
+structure = NonnegativeNeuralLyapunov(
+    dim_output;
+    δ = 1e-6
+)
+minimization_condition = DontCheckNonnegativity(check_fixed_point = true)
+
+# Define Lyapunov decrease condition
+# Damped SHO has exponential decrease at a rate of k = ζ * ω_0, so we train to certify that
+decrease_condition = ExponentialDecrease(prod(p))
+
+# Construct neural Lyapunov specification
+spec = NeuralLyapunovSpecification(
+    structure,
+    minimization_condition,
+    decrease_condition
+)
+
+# Construct PDESystem
+@named pde_system = NeuralLyapunovPDESystem(
+    dynamics,
+    lb,
+    ub,
+    spec;
+    p = p
+)
+```
+
+Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem.
+
+```@example SHO
+prob = discretize(pde_system, discretization)
+
+using Optimization, OptimizationOptimisers, OptimizationOptimJL
+
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500)
+prob = Optimization.remake(prob, u0 = res.u)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 500)
+
+net = discretization.phi
+θ = res.u.depvar
+```
+
+We can use the result of the optimization problem to build the Lyapunov candidate as a Julia function.
+
+```@example SHO
+V, V̇ = get_numerical_lyapunov_function(
+    net,
+    θ,
+    structure,
+    f,
+    fixed_point;
+    p = p
+)
+```
+
+Now let's see how we did.
+We'll evaluate both ``V`` and ``\dot{V}`` on a ``101 \times 101`` grid:
+
+```@example SHO
+Δx = (ub[1] - lb[1]) / 100
+Δv = (ub[2] - lb[2]) / 100
+xs = lb[1]:Δx:ub[1]
+vs = lb[2]:Δv:ub[2]
+states = Iterators.map(collect, Iterators.product(xs, vs))
+V_samples = vec(V(hcat(states...)))
+V̇_samples = vec(V̇(hcat(states...)))
+
+# Print statistics
+V_min, i_min = findmin(V_samples)
+state_min = collect(states)[i_min]
+V_min, state_min = if V(fixed_point) ≤ V_min
+        V(fixed_point), fixed_point
+    else
+        V_min, state_min
+    end
+
+println("V(0.,0.) = ", V(fixed_point))
+println("V ∋ [", V_min, ", ", maximum(V_samples), "]")
+println("Minimial sample of V is at ", state_min)
+println(
+    "V̇ ∋ [",
+    minimum(V̇_samples),
+    ", ",
+    max(V̇(fixed_point), maximum(V̇_samples)),
+    "]",
+)
+```
+
+At least at these validation samples, the conditions that ``\dot{V}`` be negative semi-definite and ``V`` be minimized at the origin are nearly sastisfied.
+
+```@example SHO
+using Plots
+
+p1 = plot(xs, vs, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "v");
+p1 = scatter!([0], [0], label = "Equilibrium");
+p2 = plot(
+    xs,
+    vs,
+    V̇_samples,
+    linetype = :contourf,
+    title = "V̇",
+    xlabel = "x",
+    ylabel = "v",
+);
+p2 = scatter!([0], [0], label = "Equilibrium");
+plot(p1, p2)
+```
+
+Each sublevel set of ``V`` completely contained in the plot above has been verified as a subset of the region of attraction.
\ No newline at end of file
diff --git a/docs/src/demos/policy_search.md b/docs/src/demos/policy_search.md
new file mode 100644
index 0000000..569a6c9
--- /dev/null
+++ b/docs/src/demos/policy_search.md
@@ -0,0 +1,426 @@
+# Policy Search on the Driven Inverted Pendulum
+
+In this demonstration, we'll search for a neural network policy to stabilize the upright equilibrium of the inverted pendulum.
+
+The governing differential equation for the driven pendulum is:
+```math
+\frac{d^2 \theta}{dt^2} + 2 \zeta \omega_0 \frac{d \theta}{dt} + \omega_0^2 \sin(\theta) = \tau,
+```
+where ``\theta`` is the counterclockwise angle from the downward equilibrium, ``\zeta`` and ``\omega_0`` are system parameters, and ``\tau`` is our control input&mdash;the torque.
+
+We'll jointly train a neural controller ``\tau = u \left( \theta, \frac{d\theta}{dt} \right)`` and neural Lyapunov function ``V`` which will certify the stability of the closed-loop system.
+
+## Copy-Pastable Code
+
+```julia
+using NeuralPDE, Lux, ModelingToolkit, NeuralLyapunov
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
+using Random
+
+Random.seed!(200)
+
+######################### Define dynamics and domain ##########################
+
+@parameters ζ ω_0
+defaults = Dict([ζ => 0.5, ω_0 => 1.0])
+
+@variables t θ(t) τ(t) [input = true]
+Dt = Differential(t)
+DDt = Dt^2
+
+eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ τ]
+
+@named driven_pendulum = ODESystem(
+    eqs,
+    t,
+    [θ, τ],
+    [ζ, ω_0];
+    defaults = defaults
+)
+
+bounds = [
+    θ ∈ (0, 2π),
+    Dt(θ) ∈ (-10.0, 10.0)
+]
+
+upright_equilibrium = [π, 0.0]
+
+####################### Specify neural Lyapunov problem #######################
+
+# Define neural network discretization
+# We use an input layer that is periodic with period 2π with respect to θ
+dim_state = length(bounds)
+dim_hidden = 20
+dim_phi = 3
+dim_u = 1
+dim_output = dim_phi + dim_u
+chain = [Lux.Chain(
+             PeriodicEmbedding([1], [2π]),
+             Dense(3, dim_hidden, tanh),
+             Dense(dim_hidden, dim_hidden, tanh),
+             Dense(dim_hidden, 1)
+         ) for _ in 1:dim_output]
+
+# Define neural network discretization
+strategy = QuasiRandomTraining(1_250)
+discretization = PhysicsInformedNN(chain, strategy)
+
+# Define neural Lyapunov structure
+structure = PositiveSemiDefiniteStructure(
+    dim_phi;
+    pos_def = function (state, fixed_point)
+        θ, ω = state
+        θ_eq, ω_eq = fixed_point
+        log(1.0 + (sin(θ) - sin(θ_eq))^2 + (cos(θ) - cos(θ_eq))^2 + (ω - ω_eq)^2)
+    end
+)
+structure = add_policy_search(
+    structure,
+    dim_u
+)
+minimization_condition = DontCheckNonnegativity(check_fixed_point = false)
+
+# Define Lyapunov decrease condition
+decrease_condition = AsymptoticDecrease(strict = true)
+
+# Construct neural Lyapunov specification
+spec = NeuralLyapunovSpecification(
+    structure,
+    minimization_condition,
+    decrease_condition
+)
+
+############################# Construct PDESystem #############################
+
+@named pde_system = NeuralLyapunovPDESystem(
+    driven_pendulum,
+    bounds,
+    spec;
+    fixed_point = upright_equilibrium
+)
+
+######################## Construct OptimizationProblem ########################
+
+prob = discretize(pde_system, discretization)
+
+########################## Solve OptimizationProblem ##########################
+
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 400)
+prob = Optimization.remake(prob, u0 = res.u)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
+
+###################### Get numerical numerical functions ######################
+
+net = discretization.phi
+_θ = res.u.depvar
+
+(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function(
+    driven_pendulum; simplify = true)
+p = [defaults[param] for param in p_order]
+
+V_func, V̇_func = get_numerical_lyapunov_function(
+    net,
+    _θ,
+    structure,
+    open_loop_pendulum_dynamics,
+    upright_equilibrium;
+    p = p
+)
+
+u = get_policy(net, _θ, dim_output, dim_u)
+```
+
+## Detailed description
+
+```@setup policy_search
+using Random
+
+Random.seed!(200)
+```
+
+In this example, we'll set up the dynamics using a ModelingToolkit `ODESystem`.
+Since the torque ``\tau`` is our control input, we use the `[input = true]` flag for it.
+
+Since the angle ``\theta`` is periodic with period ``2\pi``, our box domain will be one period in ``\theta`` and an interval in ``\frac{d\theta}{dt}``.
+
+```@example policy_search
+using ModelingToolkit
+
+@parameters ζ ω_0
+defaults = Dict([ζ => 0.5, ω_0 => 1.0])
+
+@variables t θ(t) τ(t) [input = true]
+Dt = Differential(t)
+DDt = Dt^2
+
+eqs = [DDt(θ) + 2ζ * ω_0 * Dt(θ) + ω_0^2 * sin(θ) ~ τ]
+
+bounds = [
+    θ ∈ (0, 2π),
+    Dt(θ) ∈ (-10.0, 10.0)
+]
+
+upright_equilibrium = [π, 0.0]
+
+@named driven_pendulum = ODESystem(
+    eqs,
+    t,
+    [θ, τ],
+    [ζ, ω_0];
+    defaults = defaults
+)
+```
+
+We'll use an architecture that's ``2\pi``-periodic in ``\theta`` so that we can train on just one period of ``\theta`` and don't need to add any periodic boundary conditions.
+To achieve that, we use `Lux.PeriodicEmbedding([1], [2pi])`, enforces `2pi`-periodicity in input number `1`.
+Additionally, we include output dimensions for both the neural Lyapunov function and the neural controller.
+
+Other than that, setting up the neural network using Lux and NeuralPDE training strategy is no different from any other physics-informed neural network problem.
+For more on that aspect, see the [NeuralPDE documentation](https://docs.sciml.ai/NeuralPDE/stable/).
+
+```@example policy_search
+using Lux
+
+# Define neural network discretization
+# We use an input layer that is periodic with period 2π with respect to θ
+dim_state = length(bounds)
+dim_hidden = 20
+dim_phi = 3
+dim_u = 1
+dim_output = dim_phi + dim_u
+chain = [Lux.Chain(
+             PeriodicEmbedding([1], [2π]),
+             Dense(3, dim_hidden, tanh),
+             Dense(dim_hidden, dim_hidden, tanh),
+             Dense(dim_hidden, 1)
+         ) for _ in 1:dim_output]
+```
+
+```@example policy_search
+using NeuralPDE
+
+# Define neural network discretization
+strategy = QuasiRandomTraining(1250)
+discretization = PhysicsInformedNN(chain, strategy)
+```
+
+We now define our Lyapunov candidate structure along with the form of the Lyapunov conditions we'll be using.
+
+The default Lyapunov candidate from [`PositiveSemiDefiniteStructure`](@ref) is:
+```math
+V(x) = \left( 1 + \lVert \phi(x) \rVert^2 \right) \log \left( 1 + \lVert x \rVert^2 \right),
+```
+which structurally enforces positive definiteness.
+We'll modify the second factor to be ``2\pi``-periodic in ``\theta``:
+
+```@example policy_search
+using NeuralLyapunov
+
+# Define neural Lyapunov structure
+structure = PositiveSemiDefiniteStructure(
+    dim_phi;
+    pos_def = function (state, fixed_point)
+        θ, ω = state
+        θ_eq, ω_eq = fixed_point
+        log(1.0 + (sin(θ) - sin(θ_eq))^2 + (cos(θ) - cos(θ_eq))^2 + (ω - ω_eq)^2)
+    end
+)
+```
+
+In addition to representing the neural Lyapunov function, our neural network must also represent the controller.
+For this, we use the [`add_policy_search`](@ref) function, which tells NeuralLyapunov to expect dynamics with a control input and to treat the last `dim_u` dimensions of the neural network as the output of our controller.
+
+```@example policy_search
+structure = add_policy_search(
+    structure,
+    dim_u
+)
+```
+
+Since our Lyapunov candidate structurally enforces positive definiteness, we use [`DontCheckNonnegativity`](@ref).
+
+```@example policy_search
+minimization_condition = DontCheckNonnegativity(check_fixed_point = false)
+
+# Define Lyapunov decrease condition
+decrease_condition = AsymptoticDecrease(strict = true)
+
+# Construct neural Lyapunov specification
+spec = NeuralLyapunovSpecification(
+    structure,
+    minimization_condition,
+    decrease_condition
+)
+
+# Construct PDESystem 
+@named pde_system = NeuralLyapunovPDESystem(
+    driven_pendulum,
+    bounds,
+    spec;
+    fixed_point = upright_equilibrium
+)
+```
+
+Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem.
+
+```@example policy_search
+prob = discretize(pde_system, discretization)
+
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
+
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 400)
+prob = Optimization.remake(prob, u0 = res.u)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
+
+net = discretization.phi
+_θ = res.u.depvar
+```
+
+We can use the result of the optimization problem to build the Lyapunov candidate as a Julia function, as well as extract our controller, using the [`get_policy`](@ref) function.
+
+```@example policy_search
+(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function(driven_pendulum; simplify = true)
+p = [defaults[param] for param in p_order]
+
+V_func, V̇_func = get_numerical_lyapunov_function(
+    net,
+    _θ,
+    structure,
+    open_loop_pendulum_dynamics,
+    upright_equilibrium;
+    p = p
+)
+
+u = get_policy(net, _θ, dim_output, dim_u)
+```
+
+Now, let's evaluate our controller.
+First, we'll get the usual summary statistics on the Lyapunov function and plot ``V``, ``\dot{V}``, and the violations of the decrease condition.
+
+```@example policy_search
+lb = [0.0, -10.0];
+ub = [2π, 10.0];
+xs = (-2π):0.1:(2π)
+ys = lb[2]:0.1:ub[2]
+states = Iterators.map(collect, Iterators.product(xs, ys))
+V_samples = vec(V_func(hcat(states...)))
+V̇_samples = vec(V̇_func(hcat(states...)))
+
+# Print statistics
+println("V(π, 0) = ", V_func(upright_equilibrium))
+println(
+    "f([π, 0], u([π, 0])) = ",
+    open_loop_pendulum_dynamics(upright_equilibrium, u(upright_equilibrium), p, 0.0)
+)
+println(
+    "V ∋ [",
+    min(V_func(upright_equilibrium),
+    minimum(V_samples)),
+    ", ",
+    maximum(V_samples),
+    "]"
+)
+println(
+    "V̇ ∋ [",
+    minimum(V̇_samples),
+    ", ",
+    max(V̇_func(upright_equilibrium), maximum(V̇_samples)),
+    "]"
+)
+```
+
+```@example policy_search
+using Plots
+
+p1 = plot(
+    xs / pi,
+    ys,
+    V_samples,
+    linetype =
+    :contourf,
+    title = "V",
+    xlabel = "θ/π",
+    ylabel = "ω",
+    c = :bone_1
+);
+p1 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0],
+    label = "Downward Equilibria", color = :red, markershape = :x);
+p1 = scatter!(
+    [-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green, markershape = :+);
+p2 = plot(
+    xs / pi,
+    ys,
+    V̇_samples,
+    linetype = :contourf,
+    title = "dV/dt",
+    xlabel = "θ/π",
+    ylabel = "ω",
+    c = :binary
+);
+p2 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0],
+    label = "Downward Equilibria", color = :red, markershape = :x);
+p2 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green,
+    markershape = :+, legend = false);
+p3 = plot(
+    xs / pi,
+    ys,
+    V̇_samples .< 0,
+    linetype = :contourf,
+    title = "dV/dt < 0",
+    xlabel = "θ/π",
+    ylabel = "ω",
+    colorbar = false,
+    linewidth = 0
+);
+p3 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0],
+    label = "Downward Equilibria", color = :green, markershape = :+);
+p3 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria",
+    color = :red, markershape = :x, legend = false);
+plot(p1, p2, p3)
+```
+
+Now, let's simulate the closed-loop dynamics to verify that the controller can get our system to the upward equilibrium.
+
+First, we'll start at the downward equilibrium:
+
+```@example policy_search
+state_order = map(st -> SymbolicUtils.iscall(st) ? operation(st) : st, state_order)
+state_syms = Symbol.(state_order)
+
+closed_loop_dynamics = ODEFunction(
+    (x, p, t) -> open_loop_pendulum_dynamics(x, u(x), p, t);
+    sys = SciMLBase.SymbolCache(state_syms, Symbol.(p_order))
+)
+
+using DifferentialEquations
+
+# Starting still at bottom ...
+downward_equilibrium = zeros(2)
+ode_prob = ODEProblem(closed_loop_dynamics, downward_equilibrium, [0.0, 120.0], p)
+sol = solve(ode_prob, Tsit5())
+plot(sol)
+```
+
+```@example policy_search
+# ...the system should make it to the top
+θ_end, ω_end = sol.u[end]
+x_end, y_end = sin(θ_end), -cos(θ_end)
+[x_end, y_end, ω_end] # Should be approximately [0.0, 1.0, 0.0]
+```
+
+Then, we'll start at a random state:
+
+```@example policy_search
+# Starting at a random point ...
+x0 = lb .+ rand(2) .* (ub .- lb)
+ode_prob = ODEProblem(closed_loop_dynamics, x0, [0.0, 150.0], p)
+sol = solve(ode_prob, Tsit5())
+plot(sol)
+```
+
+```@example policy_search
+# ...the system should make it to the top
+θ_end, ω_end = sol.u[end]
+x_end, y_end = sin(θ_end), -cos(θ_end)
+[x_end, y_end, ω_end] # Should be approximately [0.0, 1.0, 0.0]
+```
diff --git a/docs/src/demos/roa_estimation.md b/docs/src/demos/roa_estimation.md
new file mode 100644
index 0000000..2bbed30
--- /dev/null
+++ b/docs/src/demos/roa_estimation.md
@@ -0,0 +1,248 @@
+# Estimating the Region of Attraction
+
+In this demonstration, we add awareness of the region of attraction (RoA) estimation task to our training.
+
+We'll be examining the simple one-dimensional differential equation:
+```math
+\frac{dx}{dt} = - x + x^3.
+```
+This system has a fixed point at ``x = 0`` which has a RoA of ``x \in (-1, 1)``, which we will attempt to identify.
+
+We'll train in the larger domain ``x \in [-2, 2]``.
+
+## Copy-Pastable Code
+
+```julia
+using NeuralPDE, Lux, NeuralLyapunov
+using Optimization, OptimizationOptimisers, OptimizationOptimJL
+using Random
+
+Random.seed!(200)
+
+######################### Define dynamics and domain ##########################
+
+f(x, p, t) = -x .+ x .^ 3
+lb = [-2.0];
+ub = [ 2.0];
+fixed_point = [0.0];
+
+####################### Specify neural Lyapunov problem #######################
+
+# Define neural network discretization
+dim_state = length(lb)
+dim_hidden = 5
+dim_output = 2
+chain = [Lux.Chain(
+             Dense(dim_state, dim_hidden, tanh),
+             Dense(dim_hidden, dim_hidden, tanh),
+             Dense(dim_hidden, 1, use_bias = false)
+         ) for _ in 1:dim_output]
+
+# Define training strategy
+strategy = GridTraining(0.1)
+discretization = PhysicsInformedNN(chain, strategy)
+
+# Define neural Lyapunov structure
+structure = PositiveSemiDefiniteStructure(dim_output)
+minimization_condition = DontCheckNonnegativity()
+
+# Define Lyapunov decrease condition
+decrease_condition = make_RoA_aware(AsymptoticDecrease(strict = true))
+
+# Construct neural Lyapunov specification
+spec = NeuralLyapunovSpecification(
+    structure,
+    minimization_condition,
+    decrease_condition
+)
+
+############################# Construct PDESystem #############################
+
+@named pde_system = NeuralLyapunovPDESystem(
+    f,
+    lb,
+    ub,
+    spec
+)
+
+######################## Construct OptimizationProblem ########################
+
+prob = discretize(pde_system, discretization)
+
+########################## Solve OptimizationProblem ##########################
+
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300)
+prob = Optimization.remake(prob, u0 = res.u)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
+
+###################### Get numerical numerical functions ######################
+net = discretization.phi
+θ = res.u.depvar
+
+V, V̇ = get_numerical_lyapunov_function(
+    net,
+    θ,
+    structure,
+    f,
+    fixed_point
+)
+
+################################## Simulate ###################################
+states = lb[]:0.001:ub[]
+V_samples = vec(V(states'))
+V̇_samples = vec(V̇(states'))
+
+# Calculated RoA estimate
+ρ = decrease_condition.ρ
+RoA_states = states[vec(V(transpose(states))) .≤ ρ]
+RoA = (first(RoA_states), last(RoA_states))
+```
+
+## Detailed description
+
+```@setup RoA
+using Random
+
+Random.seed!(200)
+```
+
+In this example, we set up the dynamics as a Julia function and don't bother specifying the symbols for the variables (so ``x`` will be called the default `state1` in the PDESystem).
+
+```@example RoA
+f(x, p, t) = -x .+ x .^ 3
+lb = [-2.0];
+ub = [ 2.0];
+fixed_point = [0.0];
+nothing # hide
+```
+
+Setting up the neural network using Lux and NeuralPDE training strategy is no different from any other physics-informed neural network problem.
+For more on that aspect, see the [NeuralPDE documentation](https://docs.sciml.ai/NeuralPDE/stable/).
+Since we're only considering one dimension, training on a grid isn't so bad in this case.
+
+```@example RoA
+using Lux
+
+# Define neural network discretization
+dim_state = length(lb)
+dim_hidden = 5
+dim_output = 2
+chain = [Lux.Chain(
+             Dense(dim_state, dim_hidden, tanh),
+             Dense(dim_hidden, dim_hidden, tanh),
+             Dense(dim_hidden, 1, use_bias = false)
+         ) for _ in 1:dim_output]
+```
+
+```@example RoA
+using NeuralPDE
+
+# Define training strategy
+strategy = GridTraining(0.1)
+discretization = PhysicsInformedNN(chain, strategy)
+```
+
+We now define our Lyapunov candidate structure along with the form of the Lyapunov conditions we'll be using.
+
+For this example, let's use the default Lyapunov candidate from [`PositiveSemiDefiniteStructure`](@ref):
+```math
+V(x) = \left( 1 + \lVert \phi(x) \rVert^2 \right) \log \left( 1 + \lVert x \rVert^2 \right),
+```
+which structurally enforces positive definiteness.
+We therefore use [`DontCheckNonnegativity()`](@ref).
+
+We only require asymptotic decrease in this example, but we use [`make_RoA_aware`](@ref) to only penalize positive values of ``\dot{V}(x)`` when ``V(x) \le 1``.
+
+```@example RoA
+using NeuralLyapunov
+
+# Define neural Lyapunov structure
+structure = PositiveSemiDefiniteStructure(dim_output)
+minimization_condition = DontCheckNonnegativity()
+
+# Define Lyapunov decrease condition
+decrease_condition = make_RoA_aware(AsymptoticDecrease(strict = true))
+
+# Construct neural Lyapunov specification
+spec = NeuralLyapunovSpecification(
+    structure,
+    minimization_condition,
+    decrease_condition
+)
+
+# Construct PDESystem 
+@named pde_system = NeuralLyapunovPDESystem(
+    f,
+    lb,
+    ub,
+    spec
+)
+```
+
+Now, we solve the PDESystem using NeuralPDE the same way we would any PINN problem.
+
+```@example RoA
+prob = discretize(pde_system, discretization)
+
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
+
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300)
+prob = Optimization.remake(prob, u0 = res.u)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
+
+net = discretization.phi
+θ = res.u.depvar
+```
+
+We can use the result of the optimization problem to build the Lyapunov candidate as a Julia function, then sample on a finer grid than we trained on to find the estimated region of attraction.
+
+```@example RoA
+V, V̇ = get_numerical_lyapunov_function(
+    net,
+    θ,
+    structure,
+    f,
+    fixed_point
+)
+
+# Sample
+states = lb[]:0.001:ub[]
+V_samples = vec(V(states'))
+V̇_samples = vec(V̇(states'))
+
+# Calculate RoA estimate
+ρ = decrease_condition.ρ
+RoA_states = states[vec(V(transpose(states))) .≤ ρ]
+RoA = (first(RoA_states), last(RoA_states))
+
+# Print statistics
+println("V(0.,0.) = ", V(fixed_point))
+println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]")
+println(
+    "V̇ ∋ [",
+    minimum(V̇_samples),
+    ", ",
+    max(V̇(fixed_point), maximum(V̇_samples)),
+    "]",
+)
+println("True region of attraction: (-1, 1)")
+println("Estimated region of attraction: ", RoA)
+```
+
+The estimated region of attraction is within the true region of attraction.
+
+```@example RoA
+using Plots
+
+p_V = plot(states, V_samples, label = "V", xlabel = "x", linewidth=2);
+p_V = hline!([ρ], label = "V = ρ", legend = :top);
+p_V = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/);
+p_V = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green);
+
+p_V̇ = plot(states, V̇_samples, label = "dV/dt", xlabel = "x", linewidth=2);
+p_V̇ = hline!([0.0], label = "dV/dt = 0", legend = :top);
+p_V̇ = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/);
+p_V̇ = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green);
+
+plt = plot(p_V, p_V̇)
+```
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index 65267f6..9393401 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -2,13 +2,24 @@
 CurrentModule = NeuralLyapunov
 ```
 
-# NeuralLyapunov
+# NeuralLyapunov.jl
 
-Documentation for [NeuralLyapunov](https://github.com/SciML/NeuralLyapunov.jl).
+[NeuralLyapunov.jl](https://github.com/SciML/NeuralLyapunov.jl) is a library for searching for neural Lyapunov functions in Julia.
 
-```@index
+This package provides an API for setting up the search for a neural Lyapunov function. Such a search can be formulated as a partial differential inequality, and this library generates a [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) PDESystem to be solved using [NeuralPDE.jl](https://github.com/SciML/NeuralPDE.jl). Since the Lyapunov conditions can be formulated in several different ways and a neural Lyapunov function can be set up in many different forms, this library presents an extensible interface for users to choose how they wish to set up the search, with useful pre-built options for common setups.
+
+## Getting Started
+
+If this is your first time using the library, start by familiarizing yourself with the [components of a neural Lyapunov problem](man.md) in NeuralLyapunov.jl.
+Then, you can dive in with any of the following demonstrations (the [damped simple harmonic oscillator](demos/damped_SHO.md) is recommended to begin):
+
+```@contents
+Pages = Main.DEMONSTRATION_PAGES
+Depth = 1
 ```
 
-```@autodocs
-Modules = [NeuralLyapunov]
+When you begin to write your own neural Lyapunov code, especially if you hope to define your own neural Lyapunov formulation, you may find any of the following manual pages useful:
+
+```@contents
+Pages = Main.MANUAL_PAGES
 ```
diff --git a/docs/src/man.md b/docs/src/man.md
new file mode 100644
index 0000000..29de61b
--- /dev/null
+++ b/docs/src/man.md
@@ -0,0 +1,18 @@
+# Components of a Neural Lyapunov Problem
+
+For a candidate Lyapunov function $V(x)$ to certify the stability of an equilibrium point $x_0$ of the dynamical system $\frac{dx}{dt} = f(x(t))$, it must satisfy two conditions:
+1. The function $V$ must be uniquely minimized at $x_0$, and 
+2. The function $V$ must decrease along system trajectories (i.e., $V(x(t))$ decreases as long as $x(t)$ is a trajectory of the dynamical system).
+
+A neural Lyapunov function represents the candidate Lyapunov function $V$ using a neural network, sometimes modifying the output of the network slightly so as to enforce one of the above conditions.
+
+Thus, we specify our neural Lyapunov problems with three components, each answering a different question:
+1. How is $V$ structured in terms of the neural network?
+2. How is the minimization condition to be enforced?
+3. How is the decrease condition to be enforced?
+
+These three components are represented by the three fields of a [`NeuralLyapunovSpecification`](@ref) object.
+
+```@docs
+NeuralLyapunovSpecification
+```
diff --git a/docs/src/man/decrease.md b/docs/src/man/decrease.md
new file mode 100644
index 0000000..6152caa
--- /dev/null
+++ b/docs/src/man/decrease.md
@@ -0,0 +1,48 @@
+# Lyapunov Decrease Condition
+
+To represent the condition that the Lyapunov function ``V(x)`` must decrease along system trajectories, we typically introduce a new function ``\dot{V}(x) = \nabla V(x) \cdot f(x)``.
+This function represents the rate of change of ``V`` along system trajectories.
+That is to say, if ``x(t)`` is a trajectory defined by ``\frac{dx}{dt} = f(x)``, then ``\dot{V}(x(t)) = \frac{d}{dt} [ V(x(t)) ]``.
+It is then sufficient to show that ``\dot{V}(x)`` is negative away from the fixed point and zero at the fixed point, since a negative derivative means a decreasing function.
+
+Put mathematically, it is sufficient to require ``\dot{V}(x) < 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0``.
+We call such functions negative definite (around the fixed point ``x_0``).
+The weaker condition that ``\dot{V}(x) \le 0 \, \forall x \ne x_0`` and ``\dot{V}(x_0) = 0`` is negative *semi-*definiteness.
+
+The condition that ``\dot{V}(x_0) = 0`` is satisfied by the definition of ``\dot{V}`` and the fact that ``x_0`` is a fixed point, so we do not need to train for it.
+In cases where the dynamics have some dependence on the neural network (e.g., in [policy search](policy_search.md)), we should rather train directly for ``f(x_0) = 0``, since the minimization condition will also guarantee ``[\nabla V](x_0) = 0``, so ``\dot{V}(x_0) = 0``.
+
+NeuralLyapunov.jl provides the [`LyapunovDecreaseCondition`](@ref) struct for users to specify the form of the decrease condition to enforce through training.
+
+```@docs
+LyapunovDecreaseCondition
+```
+
+## Pre-defined decrease conditions
+
+```@docs
+AsymptoticDecrease
+ExponentialDecrease
+DontCheckDecrease
+```
+
+## Defining your own decrease condition
+
+```@meta
+CurrentModule = NeuralLyapunov
+```
+
+If a user wishes to define their own version of the decrease condition in a form other than
+``\texttt{rate\_metric}(V(x), \dot{V}(x)) \le - \texttt{strength}(x, x_0)``,
+they must define their own subtype of [`AbstractLyapunovDecreaseCondition`](@ref).
+
+```@docs
+AbstractLyapunovDecreaseCondition
+```
+
+When constructing the PDESystem, [`NeuralLyapunovPDESystem`](@ref) uses [`check_decrease`](@ref) to determine if it should include an equation equating the result of [`get_decrease_condition`](@ref) to zero.
+
+```@docs
+check_decrease
+get_decrease_condition
+```
diff --git a/docs/src/man/internals.md b/docs/src/man/internals.md
new file mode 100644
index 0000000..3a377ee
--- /dev/null
+++ b/docs/src/man/internals.md
@@ -0,0 +1,5 @@
+# Internals
+
+```@docs
+NeuralLyapunov.phi_to_net
+```
diff --git a/docs/src/man/local_lyapunov.md b/docs/src/man/local_lyapunov.md
new file mode 100644
index 0000000..16aa55b
--- /dev/null
+++ b/docs/src/man/local_lyapunov.md
@@ -0,0 +1,7 @@
+# Local Lyapunov analysis
+
+For comparison with neural Lyapunov methods, we also provide a function for local Lyapunov analysis by linearization.
+
+```@docs
+local_lyapunov
+```
diff --git a/docs/src/man/minimization.md b/docs/src/man/minimization.md
new file mode 100644
index 0000000..c7d74b7
--- /dev/null
+++ b/docs/src/man/minimization.md
@@ -0,0 +1,43 @@
+# Lyapunov Minimization Condition
+
+The condition that the Lyapunov function ``V(x)`` must be minimized uniquely at the fixed point ``x_0`` is often represented as a requirement for ``V(x)`` to be positive away from the fixed point and zero at the fixed point.
+Put mathematically, it is sufficient to require ``V(x) > 0 \, \forall x \ne x_0`` and ``V(x_0) = 0``.
+We call such functions positive definite (around the fixed point ``x_0``).
+The weaker condition that ``V(x) \ge 0 \, \forall x \ne x_0`` and ``V(x_0) = 0`` is positive *semi-*definiteness.
+
+NeuralLyapunov.jl provides the [`LyapunovMinimizationCondition`](@ref) struct for users to specify the form of the minimization condition to enforce through training.
+
+```@docs
+LyapunovMinimizationCondition     
+```
+
+## Pre-defined minimization conditions
+
+```@docs
+StrictlyPositiveDefinite
+PositiveSemiDefinite
+DontCheckNonnegativity
+```
+
+## Defining your own minimization condition
+
+```@meta
+CurrentModule = NeuralLyapunov
+```
+
+If a user wishes to define their own version of the minimization condition in a form other than
+``V(x) ≥ \texttt{strength}(x, x_0)``,
+they must define their own subtype of [`AbstractLyapunovMinimizationCondition`](@ref).
+
+```@docs
+AbstractLyapunovMinimizationCondition
+```
+
+When constructing the PDESystem, [`NeuralLyapunovPDESystem`](@ref) uses [`check_nonnegativity`](@ref) to determine if it should include an equation equating the result of [`get_minimization_condition`](@ref) to zero.
+It additionally uses [`check_minimal_fixed_point`](@ref) to determine if it should include the equation ``V(x_0) = 0``.
+
+```@docs
+check_nonnegativity
+check_minimal_fixed_point
+get_minimization_condition
+```
diff --git a/docs/src/man/pdesystem.md b/docs/src/man/pdesystem.md
new file mode 100644
index 0000000..8cd353b
--- /dev/null
+++ b/docs/src/man/pdesystem.md
@@ -0,0 +1,26 @@
+# Solving a Neural Lyapunov Problem
+
+NeuralLyapunov.jl represents neural Lyapunov problems as systems of partial differential equations, using the `ModelingToolkit.PDESystem` type.
+Such a `PDESystem` can then be solved using a physics-informed neural network through [NeuralPDE.jl](https://github.com/SciML/NeuralPDE.jl).
+
+Candidate Lyapunov functions will be trained within a box domain subset of the state space.
+
+```@docs
+NeuralLyapunovPDESystem
+```
+
+!!! warning
+
+    When using [`NeuralLyapunovPDESystem`](@ref), the Lyapuonv function structure, minimization and decrease conditions, and dynamics will all be symbolically traced to generate the resulting `PDESystem` equations.
+    In some cases tracing results in more efficient code, but in others it can result in inefficiencies or even errors.
+    
+    If the generated `PDESystem` is then used with NeuralPDE.jl, that library's parser will convert the equations into Julia functions representing the loss, which presents another opportunity for unexpected errors.
+
+## Extracting the numerical Lyapunov function
+
+We provide the following convenience function for generating the Lyapunov function after the parameters have been found.
+If the `PDESystem` was solved using NeuralPDE.jl and Optimization.jl, then the argument `phi` is a field of the output of `NeuralPDE.discretize` and the argument `θ` is `res.u.depvar` where `res` is the result of the optimization.
+
+```@docs
+get_numerical_lyapunov_function
+```
diff --git a/docs/src/man/policy_search.md b/docs/src/man/policy_search.md
new file mode 100644
index 0000000..2f8214b
--- /dev/null
+++ b/docs/src/man/policy_search.md
@@ -0,0 +1,13 @@
+# Policy Search and Network-Sependent Dynamics
+
+At times, we wish to model a component of the dynamics with a neural network.
+A common example is the policy search case, when the closed-loop dynamics include a neural network controller.
+In such cases, we consider the dynamics to take the form of ``\frac{dx}{dt} = f(x, u, p, t)``, where ``u`` is the control input/the contribution to the dynamics from the neural network.
+We provide the [`add_policy_search`](@ref) function to transform a [`NeuralLyapunovStructure`](@ref) to include training the neural network to represent not just the Lyapunov function, but also the relevant part of the dynamics.
+
+Similar to [`get_numerical_lyapunov_function`](@ref), we provide the [`get_policy`](@ref) convenience function to construct ``u(x)`` that can be combined with the open-loop dynamics ``f(x, u, p, t)`` to create closed loop dynamics ``f_{cl}(x, p, t) = f(x, u(x), p, t)``.
+
+```@docs
+add_policy_search
+get_policy
+```
diff --git a/docs/src/man/roa.md b/docs/src/man/roa.md
new file mode 100644
index 0000000..80c2ef0
--- /dev/null
+++ b/docs/src/man/roa.md
@@ -0,0 +1,13 @@
+# Training for Region of Attraction Identification
+
+Satisfying the [minimization](minimization.md) and [decrease](decrease.md) conditions within the training region (or any region around the fixed point, however small) is sufficient for proving local stability.
+In many cases, however, we desire an estimate of the region of attraction, rather than simply a guarantee of local stability.
+
+Any compact sublevel set wherein the minimization and decrease conditions are satisfied is an inner estimate of the region of attraction.
+Therefore, we can restrict training for those conditions to only within a predetermined sublevel set ``\{ x : V(x) \le \rho \}``.
+To do so, define a [`LyapunovDecreaseCondition`](@ref) as usual and then pass it through the [`make_RoA_aware`](@ref) function, which returns an analogous [`RoAAwareDecreaseCondition`](@ref).
+
+```@docs
+make_RoA_aware
+RoAAwareDecreaseCondition
+```
diff --git a/docs/src/man/structure.md b/docs/src/man/structure.md
new file mode 100644
index 0000000..fc071eb
--- /dev/null
+++ b/docs/src/man/structure.md
@@ -0,0 +1,57 @@
+# Structuring a Neural Lyapunov function
+
+Three simple neural Lyapunov function structures are provided, but users can always specify a custom structure using the [`NeuralLyapunovStructure`](@ref) struct.
+
+## Pre-defined structures
+
+The simplest structure is to train the neural network directly to be the Lyapunov function, which can be accomplished using an [`UnstructuredNeuralLyapunov`](@ref).
+
+```@docs
+UnstructuredNeuralLyapunov
+```
+
+The condition that the Lyapunov function ``V(x)`` must be minimized uniquely at the fixed point ``x_0`` is often represented as a requirement for ``V(x)`` to be positive away from the fixed point and zero at the fixed point.
+Put mathematically, it is sufficient to require ``V(x) > 0 \, \forall x \ne x_0`` and ``V(x_0) = 0``.
+We call such functions positive definite (around the fixed point ``x_0``).
+
+Two structures are provided which partially or fully enforce the minimization condition: [`NonnegativeNeuralLyapunov`](@ref), which structurally enforces ``V(x) \ge 0`` everywhere, and [`PositiveSemiDefiniteStructure`](@ref), which additionally enforces ``V(x_0) = 0``.
+
+```@docs
+NonnegativeNeuralLyapunov
+PositiveSemiDefiniteStructure
+```
+
+## Defining your own neural Lyapunov function structure
+
+To define a new structure for a neural Lyapunov function, one must specify the form of the Lyapunov candidate ``V`` and its time derivative along a trajectory ``\dot{V}``, as well as how to call the dynamics ``f``.
+Additionally, the dimensionality of the output of the neural network must be known in advance.
+
+```@docs
+NeuralLyapunovStructure
+```
+
+### Calling the dynamics
+
+Very generally, the dynamical system can be a system of ODEs ``\dot{x} = f(x, u, p, t)``, where ``u`` is a control input, ``p`` contains parameters, and ``f`` depends on the neural network in some way.
+To capture this variety, users must supply the function `f_call(dynamics::Function, phi::Function, state, p, t)`.
+
+The most common example is when `dynamics` takes the same form as an `ODEFunction`. 
+i.e., ``\dot{x} = \texttt{dynamics}(x, p, t)``.
+In that case, `f_call(dynamics, phi, state, p, t) = dynamics(state, p, t)`.
+
+Suppose instead, the dynamics were supplied as a function of state alone: ``\dot{x} = \texttt{dynamics}(x)``.
+Then, `f_call(dynamics, phi, state, p, t) = dynamics(state)`.
+
+Finally, suppose ``\dot{x} = \texttt{dynamics}(x, u, p, t)`` has a unidimensional control input that is being trained (via [policy search](policy_search.md)) to be the second output of the neural network.
+Then `f_call(dynamics, phi, state, p, t) = dynamics(state, phi(state)[2], p, t)`.
+
+Note that, despite the inclusion of the time variable ``t``, NeuralLyapunov.jl currently only supports time-invariant systems, so only `t = 0.0` is used.
+
+### Structuring ``V`` and ``\dot{V}``
+
+The Lyapunov candidate function ``V`` gets specified as a function `V(phi, state, fixed_point)`, where `phi` is the neural network as a function `phi(state)`.
+Note that this form allows ``V(x)`` to depend on the neural network evaluated at points other than just the input ``x``.
+
+The time derivative ``\dot{V}`` is similarly defined by a function `V̇(phi, J_phi, dynamics, state, params, t, fixed_point)`.
+The function `J_phi(state)` gives the Jacobian of the neural network `phi` at `state`.
+The function `dynamics` is as above (with parameters `params`). 
diff --git a/src/NeuralLyapunovPDESystem.jl b/src/NeuralLyapunovPDESystem.jl
index 065c01e..9a29bb6 100644
--- a/src/NeuralLyapunovPDESystem.jl
+++ b/src/NeuralLyapunovPDESystem.jl
@@ -4,35 +4,38 @@
 
 Construct a `ModelingToolkit.PDESystem` representing the specified neural Lyapunov problem.
 
-# Arguments
-- `dynamics`: the dynamical system being analyzed, represented as an `ODESystem` or the
-        function `f` such that `ẋ = f(x[, u], p, t)`; either way, the ODE should not depend
-        on time and only `t = 0.0` will be used
-- `bounds`: an array of domains, defining the training domain by bounding the states (and
-        derivatives, when applicable) of `dynamics`; only used when `dynamics isa
-        ODESystem`, otherwise use `lb` and `ub`.
-- `lb` and `ub`: the training domain will be ``[lb_1, ub_1]×[lb_2, ub_2]×...``; not used
-        when `dynamics isa ODESystem`, then use `bounds`.
-- `spec::NeuralLyapunovSpecification`: defines the Lyapunov function structure, as well as
-        the minimization and decrease conditions.
-- `fixed_point`: the equilibrium being analyzed; defaults to the origin.
-- `p`: the values of the parameters of the dynamical system being analyzed; defaults to
-        `SciMLBase.NullParameters()`; not used when `dynamics isa ODESystem`, then use the
-        default parameter values of `dynamics`.
-- `state_syms`: an array of the `Symbol` representing each state; not used when `dynamics
-        isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics isa
-        ODEFunction`, symbols stored there are used, unless overridden here; if not provided
-        here and cannot be inferred, `[:state1, :state2, ...]` will be used.
-- `parameter_syms`: an array of the `Symbol` representing each parameter; not used when
-        `dynamics isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics
-        isa ODEFunction`, symbols stored there are used, unless overridden here; if not
-        provided here and cannot be inferred, `[:param1, :param2, ...]` will be used.
-- `policy_search::Bool`: whether or not to include a loss term enforcing `fixed_point` to
-        actually be a fixed point; defaults to `false`; only used when `dynamics isa
-        Function && !(dynamics isa ODEFunction)`; when `dynamics isa ODEFunction`,
-        `policy_search` must be `false`, so should not be supplied; when `dynamics isa
-        ODESystem`, value inferred by the presence of unbound inputs.
-- `name`: the name of the constructed `PDESystem`
+# Positional Arguments
+  - `dynamics`: the dynamical system being analyzed, represented as an `ODESystem` or the
+    function `f` such that `ẋ = f(x[, u], p, t)`; either way, the ODE should not depend on
+    time and only `t = 0.0` will be used. (For an example of when `f` would have a `u`
+    argument, see [`add_policy_search`](@ref).)
+  - `bounds`: an array of domains, defining the training domain by bounding the states (and
+    derivatives, when applicable) of `dynamics`; only used when `dynamics isa
+    ODESystem`, otherwise use `lb` and `ub`.
+  - `lb` and `ub`: the training domain will be ``[lb_1, ub_1]×[lb_2, ub_2]×...``; not used
+    when `dynamics isa ODESystem`, then use `bounds`.
+  - `spec`: a [`NeuralLyapunovSpecification`](@ref) defining the Lyapunov function
+    structure, as well as the minimization and decrease conditions.
+
+# Keyword Arguments
+  - `fixed_point`: the equilibrium being analyzed; defaults to the origin.
+  - `p`: the values of the parameters of the dynamical system being analyzed; defaults to
+    `SciMLBase.NullParameters()`; not used when `dynamics isa ODESystem`, then use the
+    default parameter values of `dynamics`.
+  - `state_syms`: an array of the `Symbol` representing each state; not used when `dynamics
+    isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics isa ODEFunction`,
+    symbols stored there are used, unless overridden here; if not provided here and cannot
+    be inferred, `[:state1, :state2, ...]` will be used.
+  - `parameter_syms`: an array of the `Symbol` representing each parameter; not used when
+    `dynamics isa ODESystem`, then the symbols from `dynamics` are used; if `dynamics isa
+    ODEFunction`, symbols stored there are used, unless overridden here; if not provided
+    here and cannot be inferred, `[:param1, :param2, ...]` will be used.
+  - `policy_search::Bool`: whether or not to include a loss term enforcing `fixed_point` to
+    actually be a fixed point; defaults to `false`; only used when `dynamics isa Function &&
+    !(dynamics isa ODEFunction)`; when `dynamics isa ODEFunction`, `policy_search` should
+    not be supplied (as it must be false); when `dynamics isa ODESystem`, value inferred by
+    the presence of unbound inputs.
+  - `name`: the name of the constructed `PDESystem`
 """
 function NeuralLyapunovPDESystem(
         dynamics::Function,
diff --git a/src/conditions_specification.jl b/src/conditions_specification.jl
index 9de4922..9b62677 100644
--- a/src/conditions_specification.jl
+++ b/src/conditions_specification.jl
@@ -1,25 +1,25 @@
 """
-    NeuralLyapunovStructure
+    NeuralLyapunovStructure(V, V̇, f_call, network_dim)
 
 Specifies the structure of the neural Lyapunov function and its derivative.
 
-Allows the user to define the Lyapunov in terms of the neural network to structurally
-enforce Lyapunov conditions.
-`V(phi::Function, state, fixed_point)` takes in the neural network, the state, and the fixed
-point, and outputs the value of the Lyapunov function at `state`.
-`V̇(phi::Function, J_phi::Function, f::Function, state, params, t, fixed_point)` takes in the
-neural network, jacobian of the neural network, dynamics, state, parameters and time (for
-calling the dynamics, when relevant), and fixed point, and outputs the time derivative of
-the Lyapunov function at `state`.
-`f_call(dynamics::Function, phi::Function, state, p, t)` takes in the dynamics, the neural
-network, the state, the parameters of the dynamics, and time, and outputs the derivative of
-the state; this is useful for making closed-loop dynamics which depend on the neural
-network, such as in the policy search case
-`network_dim` is the dimension of the output of the neural network.
+Allows the user to define the Lyapunov in terms of the neural network, potentially
+structurally enforcing some Lyapunov conditions.
+
+# Fields
+  - `V(phi::Function, state, fixed_point)`: outputs the value of the Lyapunov function at
+    `state`.
+  - `V̇(phi::Function, J_phi::Function, dynamics::Function, state, params, t, fixed_point)`:
+    outputs the time derivative of the Lyapunov function at `state`.
+  - `f_call(dynamics::Function, phi::Function, state, p, t)`: outputs the derivative of the
+    state; this is useful for making closed-loop dynamics which depend on the neural
+    network, such as in the policy search case.
+  - `network_dim`: the dimension of the output of the neural network.
+
+`phi` and `J_phi` above are both functions of `state` alone.
 """
 struct NeuralLyapunovStructure
     V::Function
-    ∇V::Function
     V̇::Function
     f_call::Function
     network_dim::Integer
@@ -46,9 +46,17 @@ All concrete `AbstractLyapunovDecreaseCondition` subtypes should define the
 abstract type AbstractLyapunovDecreaseCondition end
 
 """
-    NeuralLyapunovSpecification
+    NeuralLyapunovSpecification(structure, minimzation_condition, decrease_condition)
+
+Specifies a neural Lyapunov problem.
 
-Specifies a neural Lyapunov problem
+# Fields
+  - `structure`: a [`NeuralLyapunovStructure`](@ref) specifying the relationship between the
+    neural network and the candidate Lyapunov function.
+  - `minimzation_condition`: an [`AbstractLyapunovMinimizationCondition`](@ref) specifying
+    how the minimization condition will be enforced.
+  - `decrease_condition`: an [`AbstractLyapunovDecreaseCondition`](@ref) specifying how the
+    decrease condition will be enforced.
 """
 struct NeuralLyapunovSpecification
     structure::NeuralLyapunovStructure
@@ -81,8 +89,12 @@ end
 """
     get_minimization_condition(cond::AbstractLyapunovMinimizationCondition)
 
-Return a function of `V`, `state`, and `fixed_point` that equals zero when the Lyapunov
-minimization condition is met and is greater than zero when it's violated.
+Return a function of ``V``, ``x``, and ``x_0`` that equals zero when the Lyapunov
+minimization condition is met for the Lyapunov candidate function ``V`` at the point ``x``,
+and is greater than zero if it's violated.
+
+Note that the first input, ``V``, is a function, so the minimization condition can depend on
+the value of the candidate Lyapunov function at multiple points.
 """
 function get_minimization_condition(cond::AbstractLyapunovMinimizationCondition)
     error("get_condition not implemented for AbstractLyapunovMinimizationCondition of " *
@@ -103,8 +115,11 @@ end
 """
     get_decrease_condition(cond::AbstractLyapunovDecreaseCondition)
 
-Return a function of `V`, `dVdt`, `state`, and `fixed_point` that is equal to zero
-when the Lyapunov decrease condition is met and is greater than zero when it is violated.
+Return a function of ``V``, ``V̇``, ``x``, and ``x_0`` that returns zero when the Lyapunov
+decrease condition is met and a value greater than zero when it is violated.
+
+Note that the first two inputs, ``V`` and ``V̇``, are functions, so the decrease condition
+can depend on the value of these functions at multiple points.
 """
 function get_decrease_condition(cond::AbstractLyapunovDecreaseCondition)
     error("get_condition not implemented for AbstractLyapunovDecreaseCondition of type " *
diff --git a/src/decrease_conditions.jl b/src/decrease_conditions.jl
index 684905e..6b74e83 100644
--- a/src/decrease_conditions.jl
+++ b/src/decrease_conditions.jl
@@ -1,31 +1,63 @@
 """
-    LyapunovDecreaseCondition(check_decrease, decrease, strength, rectifier)
+    LyapunovDecreaseCondition(check_decrease, rate_metric, strength, rectifier)
 
-Specifies the form of the Lyapunov conditions to be used; if `check_decrease`, training will
-enforce `decrease(V, dVdt) ≤ strength(state, fixed_point)`.
+Specifies the form of the Lyapunov decrease condition to be used.
 
-The inequality will be approximated by the equation
-    `rectifier(decrease(V, dVdt) - strength(state, fixed_point)) = 0.0`.
+# Fields
+  - `check_decrease::Bool`: whether or not to train for negativity/nonpositivity of
+    ``V̇(x)``.
+  - `rate_metric::Function`: should increase with ``V̇(x)``; used when
+    `check_decrease == true`.
+  - `strength::Function`: specifies the level of strictness for negativity training; should
+    be zero when the two inputs are equal and nonnegative otherwise; used when
+    `check_decrease == true`.
+  - `rectifier::Function`: positive when the input is positive and (approximately) zero when
+    the input is negative.
 
-If the dynamics truly have a fixed point at `fixed_point` and `dVdt` has been defined
-properly in terms of the dynamics, then `dVdt(fixed_point)` will be `0` and there is no need
-to enforce `dVdt(fixed_point) = 0`, so `check_fixed_point` defaults to `false`.
+# Training conditions
+
+If `check_decrease == true`, training will enforce:
+
+``\\texttt{rate\\_metric}(V(x), V̇(x)) ≤ -\\texttt{strength}(x, x_0).``
+
+The inequality will be approximated by the equation:
+
+``\\texttt{rectifier}(\\texttt{rate\\_metric}(V(x), V̇(x)) + \\texttt{strength}(x, x_0)) = 0.``
+
+Note that the approximate equation and inequality are identical when
+``\\texttt{rectifier}(t) = \\max(0, t)``.
+
+If the dynamics truly have a fixed point at ``x_0`` and ``V̇(x)`` is truly the rate of
+decrease of ``V(x)`` along the dynamics, then ``V̇(x_0)`` will be ``0`` and there is no need
+to train for ``V̇(x_0) = 0``.
 
 # Examples:
 
 Asymptotic decrease can be enforced by requiring
-    `dVdt ≤ -C |state - fixed_point|^2`,
+    ``V̇(x) ≤ -C \\lVert x - x_0 \\rVert^2``,
+for some positive ``C``, which corresponds to
+
+    rate_metric = (V, dVdt) -> dVdt
+    strength = (x, x0) -> C * (x - x0) ⋅ (x - x0)
+
+This can also be accomplished with [`AsymptoticDecrease`](@ref).
+
+Exponential decrease of rate ``k`` is proven by
+    ``V̇(x) ≤ - k V(x)``,
 which corresponds to
-    `decrease = (V, dVdt) -> dVdt`
-    `strength = (x, x0) -> -C * (x - x0) ⋅ (x - x0)`
 
-Exponential decrease of rate `k` is proven by `dVdt ≤ - k * V`, so corresponds to
-    `decrease = (V, dVdt) -> dVdt + k * V`
-    `strength = (x, x0) -> 0.0`
+    rate_metric = (V, dVdt) -> dVdt + k * V
+    strength = (x, x0) -> 0.0
+
+This can also be accomplished with [`ExponentialDecrease`](@ref).
+
+In either case, the rectified linear unit `rectifier = (t) -> max(zero(t), t)` exactly
+represents the inequality, but differentiable approximations of this function may be
+employed.
 """
 struct LyapunovDecreaseCondition <: AbstractLyapunovDecreaseCondition
     check_decrease::Bool
-    decrease::Function
+    rate_metric::Function
     strength::Function
     rectifier::Function
 end
@@ -37,7 +69,7 @@ end
 function get_decrease_condition(cond::LyapunovDecreaseCondition)
     if cond.check_decrease
         return (V, dVdt, x, fixed_point) -> cond.rectifier(
-            cond.decrease(V(x), dVdt(x)) - cond.strength(x, fixed_point)
+            cond.rate_metric(V(x), dVdt(x)) + cond.strength(x, fixed_point)
         )
     else
         return nothing
@@ -47,12 +79,15 @@ end
 """
     AsymptoticDecrease(; strict, C, rectifier)
 
-Construct a `LyapunovDecreaseCondition` corresponding to asymptotic decrease.
+Construct a [`LyapunovDecreaseCondition`](@ref) corresponding to asymptotic decrease.
 
-If `strict` is `false`, the condition is `dV/dt ≤ 0`, and if `strict` is `true`, the
-condition is `dV/dt ≤ - C | state - fixed_point |^2`.
+If `strict == false`, the decrease condition is
+``\\dot{V}(x) ≤ 0``,
+and if `strict == true`, the condition is
+``\\dot{V}(x) ≤ - C \\lVert x - x_0 \\rVert^2``.
 
-The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`.
+The inequality is represented by
+``\\texttt{rectifier}(\\dot{V}(x) + C \\lVert x - x_0 \\rVert^2) = 0``.
 """
 function AsymptoticDecrease(;
         strict::Bool = false,
@@ -60,7 +95,7 @@ function AsymptoticDecrease(;
         rectifier = (t) -> max(zero(t), t)
 )::LyapunovDecreaseCondition
     strength = if strict
-        (x, x0) -> -C * (x - x0) ⋅ (x - x0)
+        (x, x0) -> C * (x - x0) ⋅ (x - x0)
     else
         (x, x0) -> 0.0
     end
@@ -76,12 +111,14 @@ end
 """
     ExponentialDecrease(k; strict, C, rectifier)
 
-Construct a `LyapunovDecreaseCondition` corresponding to exponential decrease of rate `k`.
+Construct a [`LyapunovDecreaseCondition`](@ref) corresponding to exponential decrease of
+rate ``k``.
 
-If `strict` is `false`, the condition is `dV/dt ≤ -k * V`, and if `strict` is `true`, the
-condition is `dV/dt ≤ -k * V - C * ||state - fixed_point||^2`.
+If `strict == false`, the condition is ``\\dot{V}(x) ≤ -k V(x)``, and if `strict == true`,
+the condition is ``\\dot{V}(x) ≤ -k V(x) - C \\lVert x - x_0 \\rVert^2``.
 
-The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`.
+The inequality is represented by
+``\\texttt{rectifier}(\\dot{V}(x) + k V(x) + C \\lVert x - x_0 \\rVert^2) = 0``.
 """
 function ExponentialDecrease(
         k::Real;
@@ -90,7 +127,7 @@ function ExponentialDecrease(
         rectifier = (t) -> max(zero(t), t)
 )::LyapunovDecreaseCondition
     strength = if strict
-        (x, x0) -> -C * (x - x0) ⋅ (x - x0)
+        (x, x0) -> C * (x - x0) ⋅ (x - x0)
     else
         (x, x0) -> 0.0
     end
@@ -106,7 +143,7 @@ end
 """
     DontCheckDecrease()
 
-Construct a `LyapunovDecreaseCondition` which represents not checking for
+Construct a [`LyapunovDecreaseCondition`](@ref) which represents not checking for
 decrease of the Lyapunov function along system trajectories. This is appropriate
 in cases when the Lyapunov decrease condition has been structurally enforced.
 """
diff --git a/src/decrease_conditions_RoA_aware.jl b/src/decrease_conditions_RoA_aware.jl
index c8ae130..321394a 100644
--- a/src/decrease_conditions_RoA_aware.jl
+++ b/src/decrease_conditions_RoA_aware.jl
@@ -1,44 +1,81 @@
 """
-    RoAAwareDecreaseCondition(check_decrease, decrease, strength, rectifier, ρ,
-                              out_of_RoA_penalty)
+    RoAAwareDecreaseCondition(check_decrease, rate_metric, strength, rectifier, ρ, out_of_RoA_penalty)
 
-Specifies the form of the Lyapunov conditions to be used, training for a region of
-attraction estimate of `{ x : V(x) ≤ ρ }`
+Specifies the form of the Lyapunov decrease condition to be used, training for a region of
+attraction estimate of ``\\{ x : V(x) ≤ ρ \\}``.
 
-If `check_decrease`, training will enforce
-`decrease(V(state), dVdt(state)) ≤ strength(state, fixed_point)` whenever `V(state) ≤ ρ`,
-and will instead apply
-`|out_of_RoA_penalty(V(state), dVdt(state), state, fixed_point, ρ)|^2` when `V(state) > ρ`.
+# Fields
+  - `check_decrease::Bool`: whether or not to train for negativity/nonpositivity of
+    ``V̇(x)``.
+  - `rate_metric::Function`: should increase with ``V̇(x)``; used when
+    `check_decrease == true`.
+  - `strength::Function`: specifies the level of strictness for negativity training; should
+    be zero when the two inputs are equal and nonnegative otherwise; used when
+    `check_decrease == true`.
+  - `rectifier::Function`: positive when the input is positive and (approximately) zero when
+    the input is negative.
+  - `ρ`: the level of the sublevel set forming the estimate of the region of attraction.
+  - `out_of_RoA_penalty::Function`: a loss function to be applied penalizing points outside
+    the sublevel set forming the region of attraction estimate.
+
+# Training conditions
+
+If `check_decrease == true`, training will enforce
+
+``\\texttt{rate\\_metric}(V(x), V̇(x)) ≤ - \\texttt{strength}(x, x_0)``
+
+whenever ``V(x) ≤ ρ``, and will instead apply a loss of
+
+``\\lvert \\texttt{out\\_of\\_RoA\\_penalty}(V(x), V̇(x), x, x_0, ρ) \\rvert^2``
+
+when ``V(x) > ρ``.
 
 The inequality will be approximated by the equation
-    `rectifier(decrease(V(state), dVdt(state)) - strength(state, fixed_point)) = 0.0`.
 
-If the dynamics truly have a fixed point at `fixed_point` and `dVdt` has been defined
-properly in terms of the dynamics, then `dVdt(fixed_point)` will be `0` and there is no need
-to enforce `dVdt(fixed_point) = 0`, so `check_fixed_point` defaults to `false`.
+``\\texttt{rectifier}(\\texttt{rate\\_metric}(V(x), V̇(x)) + \\texttt{strength}(x, x_0)) = 0``.
+
+Note that the approximate equation and inequality are identical when
+``\\texttt{rectifier}(t) = \\max(0, t)``.
+
+If the dynamics truly have a fixed point at ``x_0`` and ``V̇(x)`` is truly the rate of
+decrease of ``V(x)`` along the dynamics, then ``V̇(x_0)`` will be ``0`` and there is no need
+to train for ``V̇(x_0) = 0``.
 
 # Examples:
 
 Asymptotic decrease can be enforced by requiring
-    `dVdt ≤ -C |state - fixed_point|^2`,
+    ``V̇(x) ≤ -C \\lVert x - x_0 \\rVert^2``,
+for some positive ``C``, which corresponds to
+
+    rate_metric = (V, dVdt) -> dVdt
+    strength = (x, x0) -> C * (x - x0) ⋅ (x - x0)
+
+Exponential decrease of rate ``k`` is proven by
+    ``V̇(x) ≤ - k V(x)``,
 which corresponds to
-    `decrease = (V, dVdt) -> dVdt` and
-    `strength = (x, x0) -> -C * (x - x0) ⋅ (x - x0)`.
 
-Exponential decrease of rate `k` is proven by `dVdt ≤ - k * V`, so corresponds to
-    `decrease = (V, dVdt) -> dVdt + k * V` and
-    `strength = (x, x0) -> 0.0`.
+    rate_metric = (V, dVdt) -> dVdt + k * V
+    strength = (x, x0) -> 0.0
 
 Enforcing either condition only in the region of attraction and not penalizing any points
 outside that region would correspond to
-    `out_of_RoA_penalty = (V, dVdt, state, fixed_point, ρ) -> 0.0`,
+
+    out_of_RoA_penalty = (V, dVdt, state, fixed_point, ρ) -> 0.0
+
 whereas an example of a penalty that decays farther in state space from the fixed point is
-    `out_of_RoA_penalty = (V, dVdt, state, fixed_point, ρ) -> 1.0 / ((x - x0) ⋅ (x - x0))`.
-Note that this penalty could also depend on values of `V`, `dVdt`, and `ρ`.
+
+    out_of_RoA_penalty = (V, dVdt, x, x0, ρ) -> 1.0 / ((x - x0) ⋅ (x - x0))
+
+Note that this penalty could also depend on values of ``V`` and ``V̇`` at various points, as
+well as ``ρ``.
+
+In any of these cases, the rectified linear unit `rectifier = (t) -> max(zero(t), t)`
+exactly represents the inequality, but differentiable approximations of this function may be
+employed.
 """
 struct RoAAwareDecreaseCondition <: AbstractLyapunovDecreaseCondition
     check_decrease::Bool
-    decrease::Function
+    rate_metric::Function
     strength::Function
     rectifier::Function
     ρ::Real
@@ -53,7 +90,7 @@ function get_decrease_condition(cond::RoAAwareDecreaseCondition)
     if cond.check_decrease
         return function (V, dVdt, x, fixed_point)
             (V(x) ≤ cond.ρ) * cond.rectifier(
-                cond.decrease(V(x), dVdt(x)) - cond.strength(x, fixed_point)
+                cond.rate_metric(V(x), dVdt(x)) + cond.strength(x, fixed_point)
             ) +
             (V(x) > cond.ρ) * cond.out_of_RoA_penalty(V(x), dVdt(x), x, fixed_point,
                 cond.ρ)
@@ -67,12 +104,15 @@ end
     make_RoA_aware(cond; ρ, out_of_RoA_penalty)
 
 Add awareness of the region of attraction (RoA) estimation task to the supplied
-`LyapunovDecreaseCondition`
+[`LyapunovDecreaseCondition`](@ref).
+
+# Arguments
+  - `cond::LyapunovDecreaseCondition`: specifies the loss to be applied when ``V(x) ≤ ρ``.
+  - `ρ`: the target level such that the RoA will be ``\\{ x : V(x) ≤ ρ \\}``.
+  - `out_of_RoA_penalty::Function`: specifies the loss to be applied when ``V(x) > ρ``.
 
-`ρ` is the target level such that the RoA will be `{ x : V(x) ≤ ρ }`.
-`cond` specifies the loss applied when `V(state) ≤ ρ`, and
-`|out_of_RoA_penalty(V(state), dVdt(state), state, fixed_point, ρ)|^2` is the loss from
-`state` values such that `V(state) > ρ`.
+The loss applied to samples ``x`` such that ``V(x) > ρ`` is
+``\\lvert \\texttt{out\\_of\\_RoA\\_penalty}(V(x), V̇(x), x, x_0, ρ) \\rvert^2``.
 """
 function make_RoA_aware(
         cond::LyapunovDecreaseCondition;
@@ -81,7 +121,7 @@ function make_RoA_aware(
 )::RoAAwareDecreaseCondition
     RoAAwareDecreaseCondition(
         cond.check_decrease,
-        cond.decrease,
+        cond.rate_metric,
         cond.strength,
         cond.rectifier,
         ρ,
diff --git a/src/local_lyapunov.jl b/src/local_lyapunov.jl
index a65e571..61f48ae 100644
--- a/src/local_lyapunov.jl
+++ b/src/local_lyapunov.jl
@@ -1,5 +1,5 @@
 """
-    get_local_lyapunov(dynamics, state_dim, optimizer_factory[, jac]; fixed_point, p)
+    local_lyapunov(dynamics, state_dim, optimizer_factory[, jac]; fixed_point, p)
 
 Use semidefinite programming to derive a quadratic Lyapunov function for the linearization
 of `dynamics` around `fixed_point`. Return `(V, dV/dt)`.
@@ -19,9 +19,8 @@ If `fixed_point` is not specified, it defaults to the origin, i.e., `zeros(state
 Parameters `p` for the dynamics should be supplied when the dynamics depend on them.
 """
 function local_lyapunov(dynamics::Function, state_dim, optimizer_factory,
-                        jac::AbstractMatrix{T}; fixed_point = zeros(T, state_dim),
-                        p = SciMLBase.NullParameters()) where T <: Number
-
+        jac::AbstractMatrix{T}; fixed_point = zeros(T, state_dim),
+        p = SciMLBase.NullParameters()) where {T <: Number}
     model = JuMP.Model(optimizer_factory)
     JuMP.set_silent(model)
 
@@ -45,8 +44,8 @@ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory,
     V(states::AbstractMatrix) = mapslices(V, states, dims = [1])
 
     # Numerical gradient of Lyapunov function
-#    ∇V(state::AbstractVector) = 2 * ( Psol * (state - fixed_point) )
-#    ∇V(states::AbstractMatrix) = mapslices(∇V, states, dims = [1])
+    #    ∇V(state::AbstractVector) = 2 * ( Psol * (state - fixed_point) )
+    #    ∇V(states::AbstractMatrix) = mapslices(∇V, states, dims = [1])
 
     # Numerical time derivative of Lyapunov function
     V̇(state::AbstractVector) = 2 * dot(dynamics(state, p, 0.0), Psol, state - fixed_point)
@@ -56,8 +55,7 @@ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory,
 end
 
 function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, jac::Function;
-                        fixed_point = zeros(state_dim), p = SciMLBase.NullParameters())
-
+        fixed_point = zeros(state_dim), p = SciMLBase.NullParameters())
     A::AbstractMatrix = jac(fixed_point, p)
     return local_lyapunov(
         dynamics,
@@ -70,8 +68,7 @@ function local_lyapunov(dynamics::Function, state_dim, optimizer_factory, jac::F
 end
 
 function local_lyapunov(dynamics::Function, state_dim, optimizer_factory;
-                        fixed_point = zeros(state_dim), p = SciMLBase.NullParameters())
-
+        fixed_point = zeros(state_dim), p = SciMLBase.NullParameters())
     A::AbstractMatrix = ForwardDiff.jacobian(x -> dynamics(x, p, 0.0), fixed_point)
     return local_lyapunov(
         dynamics,
diff --git a/src/minimization_conditions.jl b/src/minimization_conditions.jl
index ad67436..9f4355e 100644
--- a/src/minimization_conditions.jl
+++ b/src/minimization_conditions.jl
@@ -1,31 +1,49 @@
 """
-    LyapunovMinimizationCondition
+    LyapunovMinimizationCondition(check_nonnegativity, strength, rectifier, check_fixed_point)
 
-Specifies the form of the Lyapunov conditions to be used.
+Specifies the form of the Lyapunov minimization condition to be used.
+
+# Fields
+  - `check_nonnegativity::Bool`: whether or not to train for positivity/nonnegativity of
+    ``V(x)``
+  - `strength::Function`: specifies the level of strictness for positivity training; should
+    be zero when the two inputs are equal and nonnegative otherwise; used when
+    `check_nonnegativity == true`
+  - `rectifier::Function`: positive when the input is positive and (approximately) zero when
+    the input is negative
+  - `check_fixed_point`: whether or not to train for ``V(x_0) = 0``.
+
+# Training conditions
+
+If `check_nonnegativity` is `true`, training will attempt to enforce:
+
+``V(x) ≥ \\texttt{strength}(x, x_0).``
+
+The inequality will be approximated by the equation:
+
+``\\texttt{rectifier}(\\texttt{strength}(x, x_0) - V(x_0)) = 0.``
+
+Note that the approximate equation and inequality are identical when
+``\\texttt{rectifier}(t) = \\max(0, t)``.
 
-If `check_nonnegativity` is `true`, training will attempt to enforce
-    `V(state) ≥ strength(state, fixed_point)`.
-The inequality will be approximated by the equation
-    `rectifier(strength(state, fixed_point) - V(state)) = 0.0`.
 If `check_fixed_point` is `true`, then training will also attempt to enforce
-    `V(fixed_point) = 0`.
+``V(x_0) = 0``.
 
 # Examples
 
-The condition that the Lyapunov function must be minimized uniquely at the fixed point can
-be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This
-could be enfored by `V(fixed_point) ≥ ||state - fixed_point||^2`, which would be
-represented, with `check_nonnegativity = true`, by
-    strength(state, fixed_point) = ||state - fixed_point||^2,
-paired with `V(fixed_point) = 0`, which can be enforced with `check_fixed_point = true`.
+When training for a strictly positive definite ``V``, an example of an appropriate `strength`
+is ``\\texttt{strength}(x, x_0) = \\lVert x - x_0 \\rVert^2``.
+This form is used in [`StrictlyPositiveDefinite`](@ref).
 
-If `V` were structured such that it is always nonnegative, then `V(fixed_point) = 0` is all
+If ``V`` were structured such that it is always nonnegative, then ``V(x_0) = 0`` is all
 that must be enforced in training for the Lyapunov function to be uniquely minimized at
-`fixed_point`. So, in that case, we would use
+``x_0``. In that case, we would use
     `check_nonnegativity = false;  check_fixed_point = true`.
+This can also be accomplished with [`DontCheckNonnegativity(true)`](@ref).
 
-In either case, `rectifier = (t) -> max(0.0, t)` exactly represents the inequality, but
-differentiable approximations of this function may be employed.
+In either case, the rectified linear unit `rectifier = (t) -> max(zero(t), t)` exactly
+represents the inequality, but differentiable approximations of this function may be
+employed.
 """
 struct LyapunovMinimizationCondition <: AbstractLyapunovMinimizationCondition
     check_nonnegativity::Bool
@@ -51,14 +69,15 @@ function get_minimization_condition(cond::LyapunovMinimizationCondition)
 end
 
 """
-    StrictlyPositiveDefinite(C; check_fixed_point, rectifier)
+    StrictlyPositiveDefinite(; C, check_fixed_point, rectifier)
 
-Construct a `LyapunovMinimizationCondition` representing
-    `V(state) ≥ C * ||state - fixed_point||^2`.
-If `check_fixed_point` is `true`, then training will also attempt to enforce
-    `V(fixed_point) = 0`.
+Construct a [`LyapunovMinimizationCondition`](@ref) representing
+    ``V(x) ≥ C \\lVert x - x_0 \\rVert^2``.
+If `check_fixed_point == true`, then training will also attempt to enforce
+    ``V(x_0) = 0``.
 
-The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`.
+The inequality is represented by
+    ``\\texttt{rectifier}(C \\lVert x - x_0 \\rVert^2 - V(x)) = 0``.
 """
 function StrictlyPositiveDefinite(;
         check_fixed_point = true,
@@ -74,15 +93,15 @@ function StrictlyPositiveDefinite(;
 end
 
 """
-    PositiveSemiDefinite(check_fixed_point)
+    PositiveSemiDefinite(; check_fixed_point)
 
-Construct a `LyapunovMinimizationCondition` representing
-    `V(state) ≥ 0`.
-If `check_fixed_point` is `true`, then training will also attempt to enforce
-    `V(fixed_point) = 0`.
+Construct a [`LyapunovMinimizationCondition`](@ref) representing
+    ``V(x) ≥ 0``.
+If `check_fixed_point == true`, then training will also attempt to enforce
+    ``V(x_0) = 0``.
 
-The inequality is represented by `a ≥ b` <==> `rectifier(b-a) = 0.0`.
-"""
+The inequality is represented by
+    ``\\texttt{rectifier}( -V(x) ) = 0``."""
 function PositiveSemiDefinite(;
         check_fixed_point = true,
         rectifier = (t) -> max(zero(t), t)
@@ -96,17 +115,19 @@ function PositiveSemiDefinite(;
 end
 
 """
-    DontCheckNonnegativity(check_fixed_point)
+    DontCheckNonnegativity(; check_fixed_point)
 
-Construct a `LyapunovMinimizationCondition` which represents not checking for nonnegativity
-of the Lyapunov function. This is appropriate in cases where this condition has been
-structurally enforced.
+Construct a [`LyapunovMinimizationCondition`](@ref) which represents not checking for
+nonnegativity of the Lyapunov function. This is appropriate in cases where this condition
+has been structurally enforced.
 
-It is still possible to check for `V(fixed_point) = 0`, even in this case, for example if
-`V` is structured to be positive for `state ≠ fixed_point`, but it is not guaranteed
-structurally that `V(fixed_point) = 0`.
+It is still possible to check for ``V(x_0) = 0``, even in this case, for example if `V` is
+structured to be positive for ``x ≠ x_0``, but it is not guaranteed structurally that
+``V(x_0) = 0`` (such as [`NonnegativeNeuralLyapunov`](@ref)). `check_fixed_point` defaults
+to `true`, since in cases where ``V(x_0) = 0`` is enforced structurally, the equation will
+reduce to `0.0 ~ 0.0`, which gets removed in most cases.
 """
-function DontCheckNonnegativity(; check_fixed_point = false)::LyapunovMinimizationCondition
+function DontCheckNonnegativity(; check_fixed_point = true)::LyapunovMinimizationCondition
     LyapunovMinimizationCondition(
         false,
         (state, fixed_point) -> 0.0,
diff --git a/src/numerical_lyapunov_functions.jl b/src/numerical_lyapunov_functions.jl
index f915d47..7ff37fe 100644
--- a/src/numerical_lyapunov_functions.jl
+++ b/src/numerical_lyapunov_functions.jl
@@ -7,28 +7,30 @@ functions representing the Lyapunov function and its time derivative: ``V(x), V
 
 These functions can operate on a state vector or columnwise on a matrix of state vectors.
 
-# Arguments
-- `phi`: the neural network, represented as `phi(x, θ)` if the neural network has a single
-        output, or a `Vector` of the same with one entry per neural network output.
-- `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the first
-        neural network output (even if there is only one), `θ[:φ2]` the parameters of the
-        second (if there are multiple), and so on.
-- `structure`: a [`NeuralLyapunovStructure`](@ref) representing the structure of the neural
-        Lyapunov function.
-- `dynamics`: the system dynamics, as a function to be used in conjunction with
-        `structure.f_call`.
-- `fixed_point`: the equilibrium point being analyzed by the Lyapunov function.
-- `p`: parameters to be passed into `dynamics`; defaults to `SciMLBase.NullParameters()`.
-- `use_V̇_structure`: when `true`, ``V̇(x)`` is calculated using `structure.V̇`; when `false`,
-        ``V̇(x)`` is calculated using `deriv` as ``\\frac{d}{dt} V(x + t f(x))`` at
-        ``t = 0``; defaults to `false`, as it is more efficient in many cases.
-- `deriv`: a function for calculating derivatives; defaults to (and expects same arguments
-        as) `ForwardDiff.derivative`; only used when `use_V̇_structure` is `false`.
-- `jac`: a function for calculating Jacobians; defaults to (and expects same arguments as)
-        `ForwardDiff.jacobian`; only used when `use_V̇_structure` is `true`.
-- `J_net`: the Jacobian of the neural network, specified as a function
-        `J_net(phi, θ, state)`; if `isnothing(J_net)` (as is the default), `J_net` will be
-        calculated using `jac`; only used when `use_V̇_structure` is `true`.
+# Positional Arguments
+  - `phi`: the neural network, represented as `phi(x, θ)` if the neural network has a single
+    output, or a `Vector` of the same with one entry per neural network output.
+  - `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the
+    first neural network output (even if there is only one), `θ[:φ2]` the parameters of the
+    second (if there are multiple), and so on.
+  - `structure`: a [`NeuralLyapunovStructure`](@ref) representing the structure of the
+    neural Lyapunov function.
+  - `dynamics`: the system dynamics, as a function to be used in conjunction with
+    `structure.f_call`.
+  - `fixed_point`: the equilibrium point being analyzed by the Lyapunov function.
+
+# Keyword Arguments
+  - `p`: parameters to be passed into `dynamics`; defaults to `SciMLBase.NullParameters()`.
+  - `use_V̇_structure`: when `true`, ``V̇(x)`` is calculated using `structure.V̇`; when `
+    false`, ``V̇(x)`` is calculated using `deriv` as ``\\frac{∂}{∂t} V(x + t f(x))`` at
+    ``t = 0``; defaults to `false`, as it is more efficient in many cases.
+  - `deriv`: a function for calculating derivatives; defaults to (and expects same arguments
+    as) `ForwardDiff.derivative`; only used when `use_V̇_structure` is `false`.
+  - `jac`: a function for calculating Jacobians; defaults to (and expects same arguments as)
+    `ForwardDiff.jacobian`; only used when `use_V̇_structure` is `true`.
+  - `J_net`: the Jacobian of the neural network, specified as a function
+    `J_net(phi, θ, state)`; if `isnothing(J_net)` (as is the default), `J_net` will be
+    calculated using `jac`; only used when `use_V̇_structure` is `true`.
 """
 function get_numerical_lyapunov_function(
         phi,
@@ -94,14 +96,13 @@ end
 Return the network as a function of state alone.
 
 # Arguments
-
-- `phi`: the neural network, represented as `phi(x, θ)` if the neural network has a single
-        output, or a `Vector` of the same with one entry per neural network output.
-- `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the first
-        neural network output (even if there is only one), `θ[:φ2]` the parameters of the
-        second (if there are multiple), and so on.
-- `idx`: the neural network outputs to include in the returned function; defaults to all and
-        only applicable when `phi isa Vector`.
+  - `phi`: the neural network, represented as `phi(x, θ)` if the neural network has a single
+    output, or a `Vector` of the same with one entry per neural network output.
+  - `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the
+    first neural network output (even if there is only one), `θ[:φ2]` the parameters of the
+    second (if there are multiple), and so on.
+  - `idx`: the neural network outputs to include in the returned function; defaults to all
+    and only applicable when `phi isa Vector`.
 """
 function phi_to_net(phi, θ)
     let _θ = θ, φ = phi
diff --git a/src/policy_search.jl b/src/policy_search.jl
index eb099b4..38a9d9c 100644
--- a/src/policy_search.jl
+++ b/src/policy_search.jl
@@ -1,25 +1,30 @@
 """
-    add_policy_search(lyapunov_structure, new_dims, control_structure)
+    add_policy_search(lyapunov_structure, new_dims; control_structure)
 
-Add dependence on the neural network to the dynamics in a `NeuralLyapunovStructure`.
+Add dependence on the neural network to the dynamics in a [`NeuralLyapunovStructure`](@ref).
 
-Add `new_dims` outputs to the neural network and feeds them through `control_structure` to
-calculate the contribution of the neural network to the dynamics.
-Use the existing `lyapunov_structure.network_dim` dimensions as in `lyapunov_structure` to
-calculate the Lyapunov function.
+# Arguments
+  - `lyapunov_structure::NeuralLyapunovStructure`: provides structure for ``V, V̇``; should
+    assume dynamics take a form of `f(x, p, t)`.
+  - `new_dims::Integer`: number of outputs of the neural network to pass into the dynamics
+    through `control_structure`.
+  - `control_structure::Function`: transforms the final `new_dims` outputs of the neural net
+    before passing them into the dynamics; defaults to passing in the neural network outputs
+    unchanged.
 
-`lyapunov_structure` should assume in its `V̇` that the dynamics take a form `f(x, p, t)`.
-The returned `NeuralLyapunovStructure` will assume instead `f(x, u, p, t)`, where `u` is the
-contribution from the neural network. Therefore, this structure cannot be used with a
-`NeuralLyapunovPDESystem` method that requires an `ODEFunction`, `ODESystem`, or
-`ODEProblem`.
+The returned `NeuralLyapunovStructure` expects dynamics of the form `f(x, u, p, t)`, where
+`u` captures the dependence of dynamics on the neural network (e.g., through a control
+input). When evaluating the dynamics, it uses `u = control_structure(phi_end(x))` where
+`phi_end` is a function that returns the final `new_dims` outputs of the neural network.
+The other `lyapunov_structure.network_dim` outputs are used for calculating ``V`` and ``V̇``,
+as specified originally by `lyapunov_structure`.
 """
 function add_policy_search(
         lyapunov_structure::NeuralLyapunovStructure,
         new_dims::Integer;
         control_structure::Function = identity
 )::NeuralLyapunovStructure
-    let V = lyapunov_structure.V, ∇V = lyapunov_structure.∇V, V̇ = lyapunov_structure.V̇,
+    let V = lyapunov_structure.V, V̇ = lyapunov_structure.V̇,
         V_dim = lyapunov_structure.network_dim, nd = new_dims, u = control_structure
 
         NeuralLyapunovStructure(
@@ -34,9 +39,6 @@ function add_policy_search(
                     V(st -> net(st)[1:V_dim, :], state, fixed_point)
                 end
             end,
-            function (net, J_net, state, fixed_point)
-                ∇V(st -> net(st)[1:V_dim], st -> J_net(st)[1:V_dim, :], state, fixed_point)
-            end,
             function (net, J_net, f, state, params, t, fixed_point)
                 V̇(st -> net(st)[1:V_dim], st -> J_net(st)[1:V_dim, :],
                     (st, p, t) -> f(st, u(net(st)[(V_dim + 1):end]), p, t), state, params,
@@ -49,16 +51,25 @@ function add_policy_search(
 end
 
 """
-    get_policy(phi, θ, network_func, dim; control_structure)
+    get_policy(phi, θ, network_dim, control_dim; control_structure)
 
-Generate a Julia function representing the control policy as a function of the state
+Generate a Julia function representing the control policy/unmodeled portion of the dynamics
+as a function of the state.
 
 The returned function can operate on a state vector or columnwise on a matrix of state
 vectors.
 
-`phi` is the neural network with parameters `θ`. The network should have `network_dim`
-outputs, the last `control_dim` of which will be passed into `control_structure` to create
-the policy output.
+# Arguments
+  - `phi`: the neural network, represented as `phi(state, θ)` if the neural network has a
+    single output, or a `Vector` of the same with one entry per neural network output.
+  - `θ`: the parameters of the neural network; `θ[:φ1]` should be the parameters of the
+    first neural network output (even if there is only one), `θ[:φ2]` the parameters of the
+    second (if there are multiple), and so on.
+  - `network_dim`: total number of neural network outputs.
+  - `control_dim`: number of neural network outputs used in the control policy.
+  - `control_structure`: transforms the final `control_dim` outputs of the neural net before
+    passing them into the dynamics; defaults to passing in the neural network outputs
+    unchanged.
 """
 function get_policy(
         phi,
diff --git a/src/structure_specification.jl b/src/structure_specification.jl
index 0b356c6..fabe589 100644
--- a/src/structure_specification.jl
+++ b/src/structure_specification.jl
@@ -1,16 +1,17 @@
 """
     UnstructuredNeuralLyapunov()
 
-Create a `NeuralLyapunovStructure` where the Lyapunov function is the neural network
+Create a [`NeuralLyapunovStructure`](@ref) where the Lyapunov function is the neural network
 evaluated at the state. This does not structurally enforce any Lyapunov conditions.
 
+Corresponds to ``V(x) = ϕ(x)``, where ``ϕ`` is the neural network.
+
 Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For
-`f(state, input, p, t)`, consider using `add_policy_search`.
+`f(state, input, p, t)`, consider using [`add_policy_search`](@ref).
 """
 function UnstructuredNeuralLyapunov()::NeuralLyapunovStructure
     NeuralLyapunovStructure(
         (net, state, fixed_point) -> net(state),
-        (net, grad_net, state, fixed_point) -> grad_net(state),
         (net, grad_net, f, state, params, t, fixed_point) -> grad_net(state) ⋅
                                                              f(state, params, t),
         (f, net, state, p, t) -> f(state, p, t),
@@ -21,25 +22,36 @@ end
 """
     NonnegativeNeuralLyapunov(network_dim, δ, pos_def; grad_pos_def, grad)
 
-Create a `NeuralLyapunovStructure` where the Lyapunov function is the L2 norm of the neural
-network output plus a constant δ times a function `pos_def`.
+Create a [`NeuralLyapunovStructure`](@ref) where the Lyapunov function is the L2 norm of the
+neural network output plus a constant δ times a function `pos_def`.
+
+Corresponds to ``V(x) = \\lVert ϕ(x) \\rVert^2 + δ \\, \\texttt{pos\\_def}(x, x_0)``, where
+``ϕ`` is the neural network and `x_0` is the equilibrium point.
 
-The condition that the Lyapunov function must be minimized uniquely at the fixed point can
-be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This
-structure ensures `V(state) ≥ 0`. Further, if `δ > 0` and
-`pos_def(fixed_point, fixed_point) = 0`, but `pos_def(state, fixed_point) > 0` when
-`state ≠ fixed_point`, this ensures that `V(state) > 0` when `state != fixed_point`. This
-does not enforce `V(fixed_point) = 0`, so that condition must included in the neural
-Lyapunov loss function.
+This structure ensures ``V(x) ≥ 0 \\, ∀ x`` when ``δ ≥ 0`` and `pos_def` is always
+nonnegative. Further, if ``δ > 0`` and `pos_def` is strictly positive definite around
+`fixed_point`, the structure ensures that ``V(x)`` is strictly positive away from
+`fixed_point`. In such cases, the minimization condition reduces to ensuring
+``V(\\texttt{fixed\\_point}) = 0``, and so [`DontCheckNonnegativity(true)`](@ref) should be
+used.
 
-`grad_pos_def(state, fixed_point)` should be the gradient of `pos_def` with respect to
-`state` at `state`. If `grad_pos_def` is not defined, it is evaluated using `grad`, which
-defaults to `ForwardDiff.gradient`.
+# Arguments
+  - `network_dim`: output dimensionality of the neural network.
+  - `δ`: weight of `pos_def`, as above.
+  - `pos_def(state, fixed_point)`: a function that is postive (semi-)definite in `state`
+    around `fixed_point`; defaults to ``log(1 + \\lVert x - x_0 \\rVert^2)``.
 
-The neural network output has dimension `network_dim`.
+# Keyword Arguments
+  - `grad_pos_def(state, fixed_point)`: the gradient of `pos_def` with respect to `state` at
+    `state`. If `isnothing(grad_pos_def)` (as is the default), the gradient of `pos_def`
+    will be evaluated using `grad`.
+  - `grad`: a function for evaluating gradients to be used when `isnothing(grad_pos_def)`;
+    defaults to, and expects the same arguments as, `ForwardDiff.gradient`.
 
 Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For
-`f(state, input, p, t)`, consider using `add_policy_search`.
+`f(state, input, p, t)`, consider using [`add_policy_search`](@ref).
+
+See also: [`DontCheckNonnegativity`](@ref)
 """
 function NonnegativeNeuralLyapunov(
         network_dim::Integer;
@@ -53,7 +65,6 @@ function NonnegativeNeuralLyapunov(
     if δ == 0.0
         NeuralLyapunovStructure(
             (net, state, fixed_point) -> net(state) ⋅ net(state),
-            (net, J_net, state, fixed_point) -> 2 * transpose(net(state)) * J_net(state),
             (net, J_net, f, state, params, t, fixed_point) -> 2 *
                                                               dot(
                 net(state), J_net(state), f(state, params, t)),
@@ -69,8 +80,6 @@ function NonnegativeNeuralLyapunov(
         NeuralLyapunovStructure(
             (net, state, fixed_point) -> net(state) ⋅ net(state) +
                                          δ * pos_def(state, fixed_point),
-            (net, J_net, state, fixed_point) -> 2 * transpose(net(state)) * J_net(state) +
-                                                δ * grad_pos_def(state, fixed_point),
             (net, J_net, f, state, params, t, fixed_point) -> 2 * dot(
                 net(state),
                 J_net(state),
@@ -85,27 +94,42 @@ end
 """
     PositiveSemiDefiniteStructure(network_dim; pos_def, non_neg, grad_pos_def, grad_non_neg, grad)
 
-Create a `NeuralLyapunovStructure` where the Lyapunov function is the product of a positive
-(semi-)definite function `pos_def` which does not depend on the network and a nonnegative
-function non_neg which does depend the network.
+Create a [`NeuralLyapunovStructure`](@ref) where the Lyapunov function is the product of a
+positive (semi-)definite function `pos_def` which does not depend on the network and a
+nonnegative function `non_neg` which does depend the network.
 
-The condition that the Lyapunov function must be minimized uniquely at the fixed point can
-be represented as `V(fixed_point) = 0`, `V(state) > 0` when `state ≠ fixed_point`. This
-structure ensures `V(state) ≥ 0`. Further, if `pos_def` is `0` only at `fixed_point` (and
-positive elsewhere) and if `non_neg` is strictly positive away from `fixed_point` (as is the
-case for the default values of `pos_def` and `non_neg`), then this structure ensures
-`V(fixed_point) = 0` and `V(state) > 0` when `state ≠ fixed_point`.
+Corresponds to ``V(x) = \\texttt{pos\\_def}(x, x_0) * \\texttt{non\\_neg}(ϕ, x, x_0)``, where
+``ϕ`` is the neural network and `x_0` is the equilibrium point.
 
-`grad_pos_def(state, fixed_point)` should be the gradient of `pos_def` with respect to
-`state` at `state`. Similarly, `grad_non_neg(net, J_net, state, fixed_point)` should be the
-gradient of `non_neg(net, state, fixed_point)` with respect to `state` at `state`. If
-`grad_pos_def` or `grad_non_neg` is not defined, it is evaluated using `grad`, which
-defaults to `ForwardDiff.gradient`.
+This structure ensures ``V(x) ≥ 0``. Further, if `pos_def` is strictly positive definite
+`fixed_point` and `non_neg` is strictly positive (as is the case for the default values of
+`pos_def` and `non_neg`), then this structure ensures ``V(x)`` is strictly positive definite
+around `fixed_point`. In such cases, the minimization condition is satisfied structurally,
+so [`DontCheckNonnegativity(false)`](@ref) should be used.
 
-The neural network output has dimension `network_dim`.
+# Arguments
+  - network_dim: output dimensionality of the neural network.
+
+# Keyword Arguments
+  - `pos_def(state, fixed_point)`: a function that is postive (semi-)definite in `state`
+    around `fixed_point`; defaults to ``log(1 + \\lVert x - x_0 \\rVert^2)``.
+  - `non_neg(net, state, fixed_point)`: a nonnegative function of the neural network; note
+    that `net` is the neural network ``ϕ``, and `net(state)` is the value of the neural
+    network at a point ``ϕ(x)``; defaults to ``1 + \\lVert ϕ(x) \\rVert^2``.
+  - `grad_pos_def(state, fixed_point)`: the gradient of `pos_def` with respect to `state` at
+    `state`. If `isnothing(grad_pos_def)` (as is the default), the gradient of `pos_def`
+    will be evaluated using `grad`.
+  - `grad_non_neg(net, J_net, state, fixed_point)`: the gradient of `non_neg` with respect
+    to `state` at `state`; `J_net` is a function outputting the Jacobian of `net` at the
+    input. If `isnothing(grad_non_neg)` (as is the default), the gradient of `non_neg` will
+    be evaluated using `grad`.
+  - `grad`: a function for evaluating gradients to be used when `isnothing(grad_pos_def) ||
+    isnothing(grad_non_neg)`; defaults to, and expects the same arguments as, `ForwardDiff.gradient`.
 
 Dynamics are assumed to be in `f(state, p, t)` form, as in an `ODEFunction`. For
-`f(state, input, p, t)`, consider using `add_policy_search`.
+`f(state, input, p, t)`, consider using [`add_policy_search`](@ref).
+
+See also: [`DontCheckNonnegativity`](@ref)
 """
 function PositiveSemiDefiniteStructure(
         network_dim::Integer;
@@ -133,10 +157,6 @@ function PositiveSemiDefiniteStructure(
     NeuralLyapunovStructure(
         (net, state, fixed_point) -> pos_def(state, fixed_point) *
                                      non_neg(net, state, fixed_point),
-        (net, J_net, state, fixed_point) -> grad_pos_def(state, fixed_point) *
-                                            non_neg(net, state, fixed_point) +
-                                            pos_def(state, fixed_point) *
-                                            grad_non_neg(net, J_net, state, fixed_point),
         (net, J_net, f, state, params, t, fixed_point) -> (f(state, params, t) ⋅
                                                            grad_pos_def(
             state, fixed_point)) *
diff --git a/test/damped_pendulum.jl b/test/damped_pendulum.jl
index df60974..a61472d 100644
--- a/test/damped_pendulum.jl
+++ b/test/damped_pendulum.jl
@@ -1,7 +1,5 @@
-using LinearAlgebra
-using NeuralPDE, Lux, ModelingToolkit
-using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt
-using NeuralLyapunov
+using NeuralPDE, Lux, ModelingToolkit, NeuralLyapunov
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
 using Random
 using Test
 
@@ -99,7 +97,7 @@ prob = discretize(pde_system, discretization)
 
 res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300)
 prob = Optimization.remake(prob, u0 = res.u)
-res = Optimization.solve(prob, BFGS(); maxiters = 300)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
 
 ###################### Get numerical numerical functions ######################
 
diff --git a/test/damped_sho.jl b/test/damped_sho.jl
index 429d810..83fb33c 100644
--- a/test/damped_sho.jl
+++ b/test/damped_sho.jl
@@ -1,7 +1,5 @@
-using LinearAlgebra
-using NeuralPDE, Lux, ModelingToolkit
-using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt
-using NeuralLyapunov
+using NeuralPDE, Lux, NeuralLyapunov
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
 using Random
 using Test
 
@@ -18,25 +16,26 @@ function f(state, p, t)
     vel = state[2]
     vcat(vel, -2ζ * ω_0 * vel - ω_0^2 * pos)
 end
-lb = [-2 * pi, -2.0];
-ub = [2 * pi, 2.0];
-p = [0.5, 1.0]
+lb = [-5.0, -2.0];
+ub = [5.0, 2.0];
+p = [0.5, 1.0];
+fixed_point = [0.0, 0.0];
 dynamics = ODEFunction(f; sys = SciMLBase.SymbolCache([:x, :v], [:ζ, :ω_0]))
 
 ####################### Specify neural Lyapunov problem #######################
 
 # Define neural network discretization
 dim_state = length(lb)
-dim_hidden = 20
-dim_output = 5
+dim_hidden = 10
+dim_output = 3
 chain = [Lux.Chain(
              Dense(dim_state, dim_hidden, tanh),
              Dense(dim_hidden, dim_hidden, tanh),
              Dense(dim_hidden, 1)
          ) for _ in 1:dim_output]
 
-# Define neural network discretization
-strategy = GridTraining(0.05)
+# Define training strategy
+strategy = QuasiRandomTraining(1000)
 discretization = PhysicsInformedNN(chain, strategy)
 
 # Define neural Lyapunov structure
@@ -74,93 +73,78 @@ sym_prob = symbolic_discretize(pde_system, discretization)
 
 ########################## Solve OptimizationProblem ##########################
 
-res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 500)
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 450)
 prob = Optimization.remake(prob, u0 = res.u)
-res = Optimization.solve(prob, BFGS(); maxiters = 500)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
 
 ###################### Get numerical numerical functions ######################
-V_func, V̇_func = get_numerical_lyapunov_function(
+V, V̇ = get_numerical_lyapunov_function(
     discretization.phi,
     res.u.depvar,
     structure,
     f,
-    zeros(2);
+    fixed_point;
     p = p,
     use_V̇_structure = true
 )
 
 ################################## Simulate ###################################
-xs, ys = [lb[i]:0.02:ub[i] for i in eachindex(lb)]
-states = Iterators.map(collect, Iterators.product(xs, ys))
-V_predict = vec(V_func(hcat(states...)))
-dVdt_predict = vec(V̇_func(hcat(states...)))
+Δx = (ub[1] - lb[1]) / 100
+Δv = (ub[2] - lb[2]) / 100
+xs = lb[1]:Δx:ub[1]
+vs = lb[2]:Δv:ub[2]
+states = Iterators.map(collect, Iterators.product(xs, vs))
+V_samples = vec(V(hcat(states...)))
+V̇_samples = vec(V̇(hcat(states...)))
 
 #################################### Tests ####################################
 
 # Network structure should enforce nonegativeness of V
-@test min(V_func([0.0, 0.0]), minimum(V_predict)) ≥ 0.0
+V_min, i_min = findmin(V_samples)
+state_min = collect(states)[i_min]
+V_min, state_min = if V(fixed_point) ≤ V_min
+    V(fixed_point), fixed_point
+else
+    V_min, state_min
+end
+@test V_min ≥ 0.0
 
-# Trained for V's minimum to be at the fixed point
-@test V_func([0.0, 0.0])≈minimum(V_predict) atol=1e-4
-@test V_func([0.0, 0.0]) < 1e-4
+# Trained for V's minimum to be near the fixed point
+@test all(abs.(state_min .- fixed_point) .≤ [Δx, Δv])
 
 # Dynamics should result in a fixed point at the origin
-@test V̇_func([0.0, 0.0]) == 0.0
+@test V̇(fixed_point) == 0.0
 
 # V̇ should be negative almost everywhere
-@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 1e-5
+@test sum(V̇_samples .> 0) / length(V̇_samples) < 1e-4
 
 #=
-# Get RoA Estimate
-data = reshape(V_predict, (length(xs), length(ys)));
-data = vcat(data[1, :], data[end, :], data[:, 1], data[:, end]);
-ρ = minimum(data)
-
 # Print statistics
-println("V(0.,0.) = ", V_func([0.0, 0.0]))
-println("V ∋ [", min(V_func([0.0, 0.0]), minimum(V_predict)), ", ", maximum(V_predict), "]")
+println("V(0.,0.) = ", V(fixed_point))
+println("V ∋ [", V_min, ", ", maximum(V_samples), "]")
+println("Minimial sample of V is at ", state_min)
 println(
     "V̇ ∋ [",
-    minimum(dVdt_predict),
+    minimum(V̇_samples),
     ", ",
-    max(V̇_func([0.0, 0.0]), maximum(dVdt_predict)),
+    max(V̇(fixed_point), maximum(V̇_samples)),
     "]",
 )
 
 # Plot results
+using Plots
 
-p1 = plot(xs, ys, V_predict, linetype = :contourf, title = "V", xlabel = "x", ylabel = "ẋ");
+p1 = plot(xs, vs, V_samples, linetype = :contourf, title = "V", xlabel = "x", ylabel = "ẋ");
 p1 = scatter!([0], [0], label = "Equilibrium");
 p2 = plot(
     xs,
-    ys,
-    dVdt_predict,
+    vs,
+    V̇_samples,
     linetype = :contourf,
     title = "dV/dt",
     xlabel = "x",
     ylabel = "ẋ",
 );
 p2 = scatter!([0], [0], label = "Equilibrium");
-p3 = plot(
-    xs,
-    ys,
-    V_predict .< ρ,
-    linetype = :contourf,
-    title = "Estimated RoA",
-    xlabel = "x",
-    ylabel = "ẋ",
-    colorbar = false,
-);
-p4 = plot(
-    xs,
-    ys,
-    dVdt_predict .< 0,
-    linetype = :contourf,
-    title = "dV/dt < 0",
-    xlabel = "x",
-    ylabel = "ẋ",
-    colorbar = false,
-);
-p4 = scatter!([0], [0], label = "Equilibrium");
-plot(p1, p2, p3, p4)
+plot(p1, p2)
 =#
diff --git a/test/inverted_pendulum.jl b/test/inverted_pendulum.jl
index ef9e274..419b050 100644
--- a/test/inverted_pendulum.jl
+++ b/test/inverted_pendulum.jl
@@ -1,7 +1,5 @@
-using LinearAlgebra
-using NeuralPDE, Lux, ModelingToolkit
-using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt
-using NeuralLyapunov
+using NeuralPDE, Lux, NeuralLyapunov
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
 using Random
 using Test
 
@@ -94,7 +92,7 @@ prob = discretize(pde_system, discretization)
 
 res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300)
 prob = Optimization.remake(prob, u0 = res.u)
-res = Optimization.solve(prob, BFGS(); maxiters = 300)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
 
 ###################### Get numerical numerical functions ######################
 
diff --git a/test/inverted_pendulum_ODESystem.jl b/test/inverted_pendulum_ODESystem.jl
index 03f02bb..de161d3 100644
--- a/test/inverted_pendulum_ODESystem.jl
+++ b/test/inverted_pendulum_ODESystem.jl
@@ -1,7 +1,5 @@
-using LinearAlgebra
-using NeuralPDE, Lux, ModelingToolkit
-using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt
-using NeuralLyapunov
+using NeuralPDE, Lux, ModelingToolkit, NeuralLyapunov
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
 using Random
 using Test
 
@@ -33,19 +31,15 @@ bounds = [
     Dt(θ) ∈ (-10.0, 10.0)
 ]
 
-(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function(
-    driven_pendulum; simplify = true)
-
 upright_equilibrium = [π, 0.0]
-p = [defaults[param] for param in p_order]
 
 ####################### Specify neural Lyapunov problem #######################
 
 # Define neural network discretization
 # We use an input layer that is periodic with period 2π with respect to θ
 dim_state = length(bounds)
-dim_hidden = 15
-dim_phi = 2
+dim_hidden = 20
+dim_phi = 3
 dim_u = 1
 dim_output = dim_phi + dim_u
 chain = [Lux.Chain(
@@ -56,7 +50,7 @@ chain = [Lux.Chain(
          ) for _ in 1:dim_output]
 
 # Define neural network discretization
-strategy = GridTraining(0.1)
+strategy = QuasiRandomTraining(1_250)
 discretization = PhysicsInformedNN(chain, strategy)
 
 # Define neural Lyapunov structure
@@ -100,22 +94,29 @@ prob = discretize(pde_system, discretization)
 
 ########################## Solve OptimizationProblem ##########################
 
-res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300)
+res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 450)
 prob = Optimization.remake(prob, u0 = res.u)
-res = Optimization.solve(prob, BFGS(); maxiters = 300)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
 
 ###################### Get numerical numerical functions ######################
 
+net = discretization.phi
+_θ = res.u.depvar
+
+(open_loop_pendulum_dynamics, _), state_order, p_order = ModelingToolkit.generate_control_function(
+    driven_pendulum; simplify = true)
+p = [defaults[param] for param in p_order]
+
 V_func, V̇_func = get_numerical_lyapunov_function(
-    discretization.phi,
-    res.u.depvar,
+    net,
+    _θ,
     structure,
     open_loop_pendulum_dynamics,
     upright_equilibrium;
     p = p
 )
 
-u = get_policy(discretization.phi, res.u.depvar, dim_output, dim_u)
+u = get_policy(net, _θ, dim_output, dim_u)
 
 ################################## Simulate ###################################
 
@@ -124,14 +125,14 @@ ub = [2π, 10.0];
 xs = (-2π):0.1:(2π)
 ys = lb[2]:0.1:ub[2]
 states = Iterators.map(collect, Iterators.product(xs, ys))
-V_predict = vec(V_func(hcat(states...)))
-dVdt_predict = vec(V̇_func(hcat(states...)))
+V_samples = vec(V_func(hcat(states...)))
+V̇_samples = vec(V̇_func(hcat(states...)))
 
 #################################### Tests ####################################
 
 # Network structure should enforce positive definiteness
 @test V_func(upright_equilibrium) == 0.0
-@test min(V_func(upright_equilibrium), minimum(V_predict)) ≥ 0.0
+@test min(V_func(upright_equilibrium), minimum(V_samples)) ≥ 0.0
 
 # Network structure should enforce periodicity in θ
 x0 = (ub .- lb) .* rand(2, 100) .+ lb
@@ -140,17 +141,17 @@ x0 = (ub .- lb) .* rand(2, 100) .+ lb
 # Training should result in a fixed point at the upright equilibrium
 @test all(isapprox.(
     open_loop_pendulum_dynamics(upright_equilibrium, u(upright_equilibrium), p, 0.0),
-    0.0; atol = 1e-8))
+    0.0; atol = 1e-3))
 @test V̇_func(upright_equilibrium) == 0.0
 
 # V̇ should be negative almost everywhere
-@test sum(dVdt_predict .> 0) / length(dVdt_predict) < 1.04e-3
+@test sum(V̇_samples .> 0) / length(V_samples) < 0.01
 
 ################################## Simulate ###################################
 
 using DifferentialEquations
 
-state_order = map(st -> istree(st) ? operation(st) : st, state_order)
+state_order = map(st -> SymbolicUtils.iscall(st) ? operation(st) : st, state_order)
 state_syms = Symbol.(state_order)
 
 closed_loop_dynamics = ODEFunction(
@@ -158,24 +159,24 @@ closed_loop_dynamics = ODEFunction(
     sys = SciMLBase.SymbolCache(state_syms, Symbol.(p_order))
 )
 
-# Starting still at bottom
+# Starting still at bottom ...
 downward_equilibrium = zeros(2)
-ode_prob = ODEProblem(closed_loop_dynamics, downward_equilibrium, [0.0, 70.0], p)
+ode_prob = ODEProblem(closed_loop_dynamics, downward_equilibrium, [0.0, 120.0], p)
 sol = solve(ode_prob, Tsit5())
 # plot(sol)
 
-# Should make it to the top
+# ...the system should make it to the top
 θ_end, ω_end = sol.u[end]
 x_end, y_end = sin(θ_end), -cos(θ_end)
 @test all(isapprox.([x_end, y_end, ω_end], [0.0, 1.0, 0.0]; atol = 1e-3))
 
-# Starting at a random point
+# Starting at a random point ...
 x0 = lb .+ rand(2) .* (ub .- lb)
-ode_prob = ODEProblem(closed_loop_dynamics, x0, [0.0, 100.0], p)
+ode_prob = ODEProblem(closed_loop_dynamics, x0, [0.0, 150.0], p)
 sol = solve(ode_prob, Tsit5())
 # plot(sol)
 
-# Should make it to the top
+# ...the system should make it to the top
 θ_end, ω_end = sol.u[end]
 x_end, y_end = sin(θ_end), -cos(θ_end)
 @test all(isapprox.([x_end, y_end, ω_end], [0.0, 1.0, 0.0]; atol = 1e-3))
@@ -190,16 +191,16 @@ println(
 println(
     "V ∋ [",
     min(V_func(upright_equilibrium),
-    minimum(V_predict)),
+    minimum(V_samples)),
     ", ",
-    maximum(V_predict),
+    maximum(V_samples),
     "]"
 )
 println(
     "V̇ ∋ [",
-    minimum(dVdt_predict),
+    minimum(V̇_samples),
     ", ",
-    max(V̇_func(upright_equilibrium), maximum(dVdt_predict)),
+    max(V̇_func(upright_equilibrium), maximum(V̇_samples)),
     "]"
 )
 
@@ -209,7 +210,7 @@ using Plots
 p1 = plot(
     xs / pi,
     ys,
-    V_predict,
+    V_samples,
     linetype =
     :contourf,
     title = "V",
@@ -218,13 +219,13 @@ p1 = plot(
     c = :bone_1
 );
 p1 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0],
-    label = "Downward Equilibria", color = :green, markershape = :+);
+    label = "Downward Equilibria", color = :red, markershape = :x);
 p1 = scatter!(
-    [-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :red, markershape = :x);
+    [-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green, markershape = :+);
 p2 = plot(
     xs / pi,
     ys,
-    dVdt_predict,
+    V̇_samples,
     linetype = :contourf,
     title = "dV/dt",
     xlabel = "θ/π",
@@ -232,13 +233,13 @@ p2 = plot(
     c = :binary
 );
 p2 = scatter!([-2 * pi, 0, 2 * pi] / pi, [0, 0, 0],
-    label = "Downward Equilibria", color = :green, markershape = :+);
-p2 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria",
-    color = :red, markershape = :x, legend = false);
+    label = "Downward Equilibria", color = :red, markershape = :x);
+p2 = scatter!([-pi, pi] / pi, [0, 0], label = "Upward Equilibria", color = :green,
+    markershape = :+, legend = false);
 p3 = plot(
     xs / pi,
     ys,
-    dVdt_predict .< 0,
+    V̇_samples .< 0,
     linetype = :contourf,
     title = "dV/dt < 0",
     xlabel = "θ/π",
diff --git a/test/local_lyapunov.jl b/test/local_lyapunov.jl
index 84e70f4..7a2c0aa 100644
--- a/test/local_lyapunov.jl
+++ b/test/local_lyapunov.jl
@@ -5,7 +5,7 @@ using NeuralLyapunov, CSDP, ForwardDiff, Test, LinearAlgebra
 f1(x, p, t) = -x
 
 # Set up Jacobian
-J1(x, p) = (-1.0*I)(3)
+J1(x, p) = (-1.0 * I)(3)
 
 # Find Lyapunov function and rescale for C = 1
 _V1, _V̇1 = local_lyapunov(f1, 3, CSDP.Optimizer, J1)
@@ -25,7 +25,7 @@ function f2(state::AbstractVector, p, t)
 end
 function f2(states::AbstractMatrix, p, t)
     ζ, ω_0 = p
-    pos, vel = states[1,:], states[2,:]
+    pos, vel = states[1, :], states[2, :]
     vcat(transpose(vel), transpose(-2ζ * ω_0 * vel - ω_0^2 * pos))
 end
 p2 = [3.2, 5.1]
@@ -37,9 +37,9 @@ dV2dt = (state) -> ForwardDiff.derivative((δt) -> V2(state + δt * f2(state, p2
 # Test V̇2 = d/dt V2
 xs2 = -1:0.03:1
 states2 = hcat(Iterators.map(collect, Iterators.product(xs2, xs2))...)
-errs2 = abs.( V̇2(states2) .- dV2dt(states2))
+errs2 = abs.(V̇2(states2) .- dV2dt(states2))
 @test all(errs2 .< 1e-10)
 
 # Test V̇2 ≺ 0
-@test all( V̇2(states2) .< 0)
+@test all(V̇2(states2) .< 0)
 @test V̇2(zeros(2)) == 0.0
diff --git a/test/roa_estimation.jl b/test/roa_estimation.jl
index c3d2418..e8ab570 100644
--- a/test/roa_estimation.jl
+++ b/test/roa_estimation.jl
@@ -1,7 +1,5 @@
-using LinearAlgebra
-using NeuralPDE, Lux, ModelingToolkit
-using Optimization, OptimizationOptimisers, OptimizationOptimJL, NLopt
-using NeuralLyapunov
+using NeuralPDE, Lux, NeuralLyapunov
+import Optimization, OptimizationOptimisers, OptimizationOptimJL
 using Random
 using Test
 
@@ -12,8 +10,9 @@ println("Region of Attraction Estimation")
 ######################### Define dynamics and domain ##########################
 
 f(x, p, t) = -x .+ x .^ 3
-lb = [-2];
-ub = [2];
+lb = [-2.0];
+ub = [2.0];
+fixed_point = [0.0];
 
 ####################### Specify neural Lyapunov problem #######################
 
@@ -27,7 +26,7 @@ chain = [Lux.Chain(
              Dense(dim_hidden, 1, use_bias = false)
          ) for _ in 1:dim_output]
 
-# Define neural network discretization
+# Define training strategy
 strategy = GridTraining(0.1)
 discretization = PhysicsInformedNN(chain, strategy)
 
@@ -63,51 +62,51 @@ sym_prob = symbolic_discretize(pde_system, discretization)
 
 res = Optimization.solve(prob, OptimizationOptimisers.Adam(); maxiters = 300)
 prob = Optimization.remake(prob, u0 = res.u)
-res = Optimization.solve(prob, BFGS(); maxiters = 300)
+res = Optimization.solve(prob, OptimizationOptimJL.BFGS(); maxiters = 300)
 
 ###################### Get numerical numerical functions ######################
-V_func, V̇_func = get_numerical_lyapunov_function(
+V, V̇ = get_numerical_lyapunov_function(
     discretization.phi,
     res.u.depvar,
     structure,
     f,
-    zeros(length(lb))
+    fixed_point
 )
 
 ################################## Simulate ###################################
-states = first(lb):0.001:first(ub)
-V_predict = vec(V_func(states'))
-dVdt_predict = vec(V̇_func(states'))
+states = lb[]:0.001:ub[]
+V_samples = vec(V(states'))
+V̇_samples = vec(V̇(states'))
 
 # Calculated RoA estimate
 ρ = decrease_condition.ρ
-RoA_states = states[vec(V_func(transpose(states))) .≤ ρ]
+RoA_states = states[vec(V(transpose(states))) .≤ ρ]
 RoA = (first(RoA_states), last(RoA_states))
 
 #################################### Tests ####################################
 
 # Network structure should enforce positive definiteness of V
-@test min(V_func([0.0]), minimum(V_predict)) ≥ 0.0
-@test V_func([0.0]) == 0.0
+@test min(V(fixed_point), minimum(V_samples)) ≥ 0.0
+@test V(fixed_point) == 0.0
 
 # Dynamics should result in a fixed point at the origin
-@test V̇_func([0.0]) == 0.0
+@test V̇(fixed_point) == 0.0
 
 # V̇ should be negative everywhere in the region of attraction except the fixed point
-@test all(V̇_func(transpose(RoA_states[RoA_states .!= 0.0])) .< 0)
+@test all(V̇(transpose(RoA_states[RoA_states .!= fixed_point[]])) .< 0)
 
 # The estimated region of attraction should be a subset of the real region of attraction
 @test first(RoA) ≥ -1.0 && last(RoA) ≤ 1.0
 
 #=
 # Print statistics
-println("V(0.,0.) = ", V_func([0.0]))
-println("V ∋ [", min(V_func([0.0]), minimum(V_predict)), ", ", maximum(V_predict), "]")
+println("V(0.,0.) = ", V(fixed_point))
+println("V ∋ [", min(V(fixed_point), minimum(V_samples)), ", ", maximum(V_samples), "]")
 println(
     "V̇ ∋ [",
-    minimum(dVdt_predict),
+    minimum(V̇_samples),
     ", ",
-    max(V̇_func([0.0]), maximum(dVdt_predict)),
+    max(V̇(fixed_point), maximum(V̇_samples)),
     "]",
 )
 println("True region of attraction: (-1, 1)")
@@ -116,12 +115,12 @@ println("Estimated region of attraction: ", RoA)
 # Plot results
 using Plots
 
-p_V = plot(states, V_predict, label = "V", xlabel = "x", linewidth=2);
+p_V = plot(states, V_samples, label = "V", xlabel = "x", linewidth=2);
 p_V = hline!([ρ], label = "V = ρ", legend = :top);
 p_V = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/);
 p_V = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green);
 
-p_V̇ = plot(states, dVdt_predict, label = "dV/dt", xlabel = "x", linewidth=2);
+p_V̇ = plot(states, V̇_samples, label = "dV/dt", xlabel = "x", linewidth=2);
 p_V̇ = hline!([0.0], label = "dV/dt = 0", legend = :top);
 p_V̇ = vspan!(collect(RoA); label = "Estimated Region of Attraction", color = :gray, fillstyle = :/);
 p_V̇ = vspan!([-1, 1]; label = "True Region of Attraction", opacity = 0.2, color = :green);