diff --git a/Manifest.toml b/Manifest.toml
index 13694c4..0bc4626 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -14,6 +14,18 @@ version = "2.3.0"
 [[ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 
+[[Arpack]]
+deps = ["Arpack_jll", "Libdl", "LinearAlgebra", "Logging"]
+git-tree-sha1 = "91ca22c4b8437da89b030f08d71db55a379ce958"
+uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+version = "0.5.3"
+
+[[Arpack_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg"]
+git-tree-sha1 = "5ba6c757e8feccf03a1554dfaf3e26b3cfc7fd5e"
+uuid = "68821587-b530-5797-8361-c406ea357684"
+version = "3.5.1+1"
+
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
@@ -410,7 +422,7 @@ uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
 version = "2.36.0+0"
 
 [[LinearAlgebra]]
-deps = ["Libdl"]
+deps = ["Libdl", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[LogExpFunctions]]
@@ -473,6 +485,10 @@ git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
 uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
 version = "1.3.5+1"
 
+[[OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
 [[OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
@@ -561,7 +577,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[Random]]
-deps = ["Serialization"]
+deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[RecipesBase]]
@@ -880,6 +896,10 @@ git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
 uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
 version = "0.15.1+0"
 
+[[libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
 [[libfdk_aac_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
diff --git a/Project.toml b/Project.toml
index 6602882..50b71db 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,8 +4,10 @@ authors = ["Alex Quinlan", "Henry Jacobson", "Pia Callmer", "Nicola Sabbadini"]
 version = "0.1.0"
 
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 CubicHermiteSpline = "2338ea46-7d83-48cf-b44e-cf8e0b5b9cc1"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
diff --git a/src/convergence_test.jl b/src/convergence_test.jl
new file mode 100644
index 0000000..3f80303
--- /dev/null
+++ b/src/convergence_test.jl
@@ -0,0 +1,47 @@
+include("Beam1D.jl")
+include("Diagnostics.jl")
+include("analytic_testcases.jl")
+import Plots
+import LinearAlgebra
+
+pars = (mu=x->1, EI=x->1, q=x->-(x<0.5))
+
+BCs  = Dict((0,'H')=>1,
+            (0,'G')=>2,
+            (1,'M')=>3,
+            (1,'Q')=>4)
+
+
+resolutions = [10, 20, 50, 100, 200, 400, 800]
+
+best_resolution = last(resolutions)*10
+grid = collect(LinRange(0, 1, best_resolution))
+sys   = Beam1D.System(Beam1D.Problem(pars,BCs,grid))
+best = Beam1D.solve_st(sys)
+
+
+errors = []
+for resolution in resolutions
+  grid = collect(LinRange(0,1,resolution))
+  sys  = Beam1D.System(Beam1D.Problem(pars,BCs,grid))
+  sol = Beam1D.solve_st(sys)
+  error = LinearAlgebra.norm(best(grid) - sol(grid)) / resolution
+  append!(errors, [error])
+end
+#plot
+plt = Plots.plot(resolutions, errors, xaxis=:log, yaxis=:log, seriestype=:scatter, label = "Convergence", legend=:bottomleft)
+Plots.ylabel!("error")
+Plots.xlabel!("resolution")
+
+ref1 = errors[1]*resolutions[1]*resolutions.^-1
+ref1 = reshape(ref1, length(errors), 1)
+ref2 = errors[1]*resolutions[1]^2*resolutions.^-2
+ref2 = reshape(ref2, length(errors), 1)
+
+Plots.plot!(resolutions, ref1, xaxis=:log, yaxis=:log, label = "Reference 1st Order")
+Plots.plot!(resolutions, ref2, xaxis=:log, yaxis=:log, label = "Reference 2nd Order")
+
+println("Plotting convergence. Press enter to continue.")
+
+Plots.gui(plt)
+readline()