diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 89dee632f4..16e9ab6548 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -3,7 +3,7 @@ $(DocStringExtensions.README)
 """
 module ModelingToolkit
 using PrecompileTools, Reexport
-@recompile_invalidations begin 
+@recompile_invalidations begin
     using DocStringExtensions
     using Compat
     using AbstractTrees
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index d99f5f69fa..39875d4cb4 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -196,6 +196,7 @@ function complete(sys::AbstractSystem)
 end
 
 for prop in [:eqs
+    :energies
     :tag
     :noiseeqs
     :iv
diff --git a/src/systems/pde/pdesystem.jl b/src/systems/pde/pdesystem.jl
index 048dbf304e..2a9dd0a410 100644
--- a/src/systems/pde/pdesystem.jl
+++ b/src/systems/pde/pdesystem.jl
@@ -37,6 +37,8 @@ domains = [t ∈ (0.0,1.0),
 struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
     "The equations which define the PDE"
     eqs::Any
+    "The energy functions"
+    energies::Any
     "The boundary conditions"
     bcs::Any
     "The domain for the independent variables."
@@ -88,6 +90,7 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
     @add_kwonly function PDESystem(eqs, bcs, domain, ivs, dvs,
         ps = SciMLBase.NullParameters();
         defaults = Dict(),
+        energies = [],
         systems = [],
         connector_type = nothing,
         metadata = nothing,
@@ -101,6 +104,7 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
         end
 
         eqs = eqs isa Vector ? eqs : [eqs]
+        energies = energies isa Vector ? energies : [energies]
 
         if !isnothing(analytic)
             analytic = analytic isa Vector ? analytic : [analytic]
@@ -123,7 +127,8 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
             analytic_func = analytic_func isa Dict ? analytic_func : analytic_func |> Dict
         end
 
-        new(eqs, bcs, domain, ivs, dvs, ps, defaults, connector_type, systems, analytic,
+        new(eqs, energies, bcs, domain, ivs, dvs, ps, defaults, connector_type, systems,
+            analytic,
             analytic_func, name, metadata, gui_metadata)
     end
 end
@@ -149,6 +154,7 @@ function Base.show(io::IO, ::MIME"text/plain", sys::PDESystem)
     println(io, summary(sys))
     println(io, "Equations: ", get_eqs(sys))
     println(io, "Boundary Conditions: ", get_bcs(sys))
+    println(io, "Energies: ", get_energies(sys))
     println(io, "Domain: ", get_domain(sys))
     println(io, "Dependent Variables: ", get_dvs(sys))
     println(io, "Independent Variables: ", get_ivs(sys))
diff --git a/test/pde.jl b/test/pde.jl
index 720b4541d2..1f8eecd918 100644
--- a/test/pde.jl
+++ b/test/pde.jl
@@ -5,10 +5,12 @@ using ModelingToolkit, DiffEqBase, LinearAlgebra, Test
 @constants h = 1
 @variables u(..)
 Dt = Differential(t)
-Dxx = Differential(x)^2
+Dx = Differential(x)
+Dxx = (Dx)^2
 eq = Dt(u(t, x)) ~ h * Dxx(u(t, x))
 bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
     u(t, 0) ~ 0, u(t, 1) ~ 0]
+energies = [(Dx(u(t, x)))^2]
 
 domains = [t ∈ (0.0, 1.0),
     x ∈ (0.0, 1.0)]
@@ -27,3 +29,29 @@ dt = 0:0.1:1
 # Test generated analytic_func
 @test all(pdesys.analytic_func[u(t, x)]([2], disct, discx) ≈
           analytic_function([2], disct, discx) for disct in dt, discx in dx)
+
+# Test with energies
+@named pdesys = PDESystem(eq,
+    bcs,
+    domains,
+    [t, x],
+    [u],
+    [h => 1],
+    energies = energies,
+    analytic = analytic)
+@show pdesys
+
+@test all(isequal.(independent_variables(pdesys), [t, x]))
+
+# Test with explicit empty energies
+@named pdesys = PDESystem(eq,
+    bcs,
+    domains,
+    [t, x],
+    [u],
+    [h => 1],
+    energies = [],
+    analytic = analytic)
+@show pdesys
+
+@test all(isequal.(independent_variables(pdesys), [t, x]))