diff --git a/examples/Project.toml b/examples/Project.toml
index 3671d5c..64aea47 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec"
diff --git a/examples/eigenvalues.jl b/examples/eigenvalues.jl
new file mode 100644
index 0000000..65b46c0
--- /dev/null
+++ b/examples/eigenvalues.jl
@@ -0,0 +1,53 @@
+using NewtonKrylov
+using Arpack
+
+using SparseArrays, LinearAlgebra
+using CairoMakie
+
+# ## 1D Bratu equations
+# $y′′ + λ * exp(y) = 0$
+
+# F(u) = 0
+
+function bratu!(res, y, Δx, λ)
+    N = length(y)
+    for i in 1:N
+        y_l = i == 1 ? zero(eltype(y)) : y[i - 1]
+        y_r = i == N ? zero(eltype(y)) : y[i + 1]
+        y′′ = (y_r - 2y[i] + y_l) / Δx^2
+
+        res[i] = y′′ + λ * exp(y[i]) # = 0
+    end
+    return nothing
+end
+
+function bratu(y, dx, λ)
+    res = similar(y)
+    bratu!(res, y, dx, λ)
+    return res
+end
+
+# ## Choice of parameters
+const N = 500
+const λ = 3.51382
+const dx = 1 / (N + 1) # Grid-spacing
+
+# ### Domain and Inital condition
+x = LinRange(0.0 + dx, 1.0 - dx, N)
+u₀ = sin.(x .* π)
+
+JOp = NewtonKrylov.JacobianOperator(
+    (res, u) -> bratu!(res, u, dx, λ),
+    similar(u₀),
+    copy(u₀);
+    symmetric = true
+)
+
+@assert issymmetric(collect(JOp))
+
+l, ϕ = eigs(JOp; nev = 10)
+l2, ϕ2 = eigs(collect(JOp); nev = 10)
+
+
+@time  eigs(collect(JOp); nev = 300)
+@time  eigs(JOp; nev = 300)
\ No newline at end of file
diff --git a/src/NewtonKrylov.jl b/src/NewtonKrylov.jl
index db993aa..68cfd5b 100644
--- a/src/NewtonKrylov.jl
+++ b/src/NewtonKrylov.jl
@@ -24,16 +24,19 @@ struct JacobianOperator{F, A}
     f_cache::F
     res::A
     u::A
-    function JacobianOperator(f::F, res, u) where {F}
+    symmetric::Union{Bool,Missing}
+    function JacobianOperator(f::F, res, u; symmetric=missing) where {F}
         f_cache = Enzyme.make_zero(f)
-        return new{F, typeof(u)}(f, f_cache, res, u)
+        return new{F, typeof(u)}(f, f_cache, res, u, symmetric, hermetian)
     end
 end
 
 Base.size(J::JacobianOperator) = (length(J.res), length(J.u))
 Base.eltype(J::JacobianOperator) = eltype(J.u)
 
-function mul!(out, J::JacobianOperator, v)
+LinearAlgebra.issymmetric(J::JacobianOperator) = J.symmetric
+
+function mul!(out::A, J::JacobianOperator{F,A}, v::A) where {F, A}
     # Enzyme.make_zero!(J.f_cache)
     f_cache = Enzyme.make_zero(J.f) # Stop gap until we can zero out mutable values
     autodiff(
@@ -45,6 +48,14 @@ function mul!(out, J::JacobianOperator, v)
     return nothing
 end
 
+function mul!(out, J::JacobianOperator, v)
+    _v = similar(J.u)
+    _v .= v
+    _out = similar(J.res)
+    mul!(_out, J, _v)
+    out .= _out
+end
+
 LinearAlgebra.adjoint(J::JacobianOperator) = Adjoint(J)
 LinearAlgebra.transpose(J::JacobianOperator) = Transpose(J)