diff --git a/Project.toml b/Project.toml
index 488817c..21834e0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.2.1"
 CCBlade = "e1828068-15df-11e9-03e4-ef195ea46fa4"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 FLOWMath = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -20,4 +21,4 @@ FLOWMath = "0.4"
 StaticArrays = "0.12, 1.0"
 VSPGeom = "0.6"
 WriteVTK = "1.8"
-julia = "1.5"
+julia = "1.10"
diff --git a/docs/make.jl b/docs/make.jl
index 3a2c684..892aa90 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,6 +1,7 @@
 using Documenter, VortexLattice
 
 makedocs(;
+    format = Documenter.HTML(;size_threshold = nothing, size_threshold_warn = 102000),
     modules = [VortexLattice],
     pages = [
         "Home" => "index.md",
diff --git a/docs/src/advanced.md b/docs/src/advanced.md
index 6198832..55ba4aa 100644
--- a/docs/src/advanced.md
+++ b/docs/src/advanced.md
@@ -20,6 +20,8 @@ system = System(grids; sections)
 nonlinear_analysis!(system; max_iter=1, tol=1E-6, damping=0.01, print_iters=false)
 ```
 
+In the above example, leaving sections[2] empty implies that surface[2] does not have airfoils and does not need any nonlinear analysis.
+
 **nonlinear_analysis!** populates the sections with angles of attack, coefficients of lift and drag, and force-per-unit length normalized by the density (defaulted to 1) and the reference velocity.
 
 ### Warnings
diff --git a/docs/src/library.md b/docs/src/library.md
index 4f8dd2d..d834f99 100644
--- a/docs/src/library.md
+++ b/docs/src/library.md
@@ -87,6 +87,12 @@ stability_derivatives
 write_vtk
 ```
 
+### Saving and Loading
+```@docs
+save_system_to_bson
+load_system_from_bson
+```
+
 ## Private API
 
 ### Geometry
diff --git a/src/VortexLattice.jl b/src/VortexLattice.jl
index 54298f7..5f9692a 100644
--- a/src/VortexLattice.jl
+++ b/src/VortexLattice.jl
@@ -7,6 +7,7 @@ using WriteVTK
 using VSPGeom
 using CCBlade
 using DelimitedFiles
+using JLD2
 
 # value for dimensionalizing, included just for clarity in the algorithms
 const RHO = 1.0
@@ -47,6 +48,7 @@ include("circulation.jl")
 include("system.jl")
 export System
 export PanelProperties, get_surface_properties
+export save_system_to_bson, load_system_from_bson
 
 include("analyses.jl")
 export steady_analysis, steady_analysis!
diff --git a/src/nonlinear.jl b/src/nonlinear.jl
index ef14f17..7bde371 100644
--- a/src/nonlinear.jl
+++ b/src/nonlinear.jl
@@ -76,7 +76,8 @@ end
 function redefine_gamma_index!(sections, ns, nc)
     gamma_start = 1
     for i in eachindex(sections)
-        if !isassigned(sections[i],1)
+        # if !isassigned(sections[i],1)
+        if isempty(sections[i])
             gamma_start += ns[i] * nc[i]
             continue
         end
@@ -137,7 +138,8 @@ function nonlinear_analysis!(system, ref, fs; max_iter=1, tol=1E-6, damping=0.01
     for i in eachindex(system.sections)
         surface_properties = system.properties[i]
         for j in eachindex(system.sections[i])
-            !isassigned(system.sections[i],1) && continue
+            # !isassigned(system.sections[i],1) && continue
+            isempty(system.sections[i]) && continue # Ensure sections are assigned, if not do not perform nonlinear analysis on this surface
             section = system.sections[i][j]
             Γavg = 0.0
             for k in eachindex(section.panels)
@@ -178,9 +180,8 @@ function _nonlinear_analysis!(system, r, damping, tol, vel, vx, vy, vz, x_c)
     vz .= 0.0
     x_c .= 0.0
     for i in eachindex(system.surfaces)
+        isempty(system.sections[i]) && continue # Ensure sections are assigned, if not do not perform nonlinear analysis on this surface
         sections = system.sections[i]
-        !isassigned(sections,1) && continue # Ensure sections are assigned, if not do not perform nonlinear analysis on this surface
-
         surface = system.surfaces[i]
         nc = size(surface, 1)
         properties = system.properties[i]
@@ -259,11 +260,11 @@ end
 
 function update_section_forces!(system, vel, vx, vy, vz, x_c)
     for i in eachindex(system.surfaces)
+        isempty(system.sections[i]) && continue # Ensure sections are assigned, if not do not perform nonlinear analysis on this surface
         surface = system.surfaces[i]
         nc = size(surface, 1)
         properties = system.properties[i]
         sections = system.sections[i]
-        !isassigned(sections,1) && continue
         vx_view = view(vx,1:nc)
         vy_view = view(vy,1:nc)
         vz_view = view(vz,1:nc)
diff --git a/src/system.jl b/src/system.jl
index 38e7623..7d2fd99 100644
--- a/src/system.jl
+++ b/src/system.jl
@@ -221,7 +221,7 @@ function System(TF::Type, nc, ns; nw = zero(nc), grids = nothing, ratios = nothi
     end
 
     if isnothing(sections)
-        sections = [Vector{SectionProperties{TF}}(undef, ns[i]) for i = 1:nsurf]
+        sections = [Vector{SectionProperties{TF}}() for i = 1:nsurf]
     else
         redefine_gamma_index!(sections, ns, nc)
     end
@@ -274,3 +274,32 @@ of panel properties (see [`PanelProperties`](@ref)) of shape (nc, ns) where `nc`
 is the number of chordwise panels and `ns` is the number of spanwise panels
 """
 get_surface_properties(system) = system.properties
+
+
+"""
+    save_system_to_bson(system, filename)
+Save the system to a BSON file.
+# Arguments:
+ - `system`: The system to save
+ - `filename`: The name of the file to save the system to
+# Returns:
+ - `nothing`: The function does not return anything, it saves the system to a file
+"""
+function save_system_to_bson(system::System, filename::AbstractString)
+    JLD2.save_object(filename, system)
+    return nothing
+end
+
+
+"""
+    load_system_from_bson(filename::AbstractString)
+Load a system from a BSON file.
+# Arguments:
+ - `filename`: The name of the file to load the system from
+# Returns:
+ - `system`: The system loaded from the file
+"""
+function load_system_from_bson(filename::AbstractString)
+    system = JLD2.load_object(filename)
+    return system
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index d6d8e4a..49db73d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1443,3 +1443,56 @@ end
     @test isapprox(CF, CF_true, atol=1e-5)
     @test isapprox(CM, CM_true, atol=1e-5)
 end
+
+@testset "save/load system" begin
+    xle = [0.0, 0.4]
+    yle = [0.0, 7.5]
+    zle = [0.0, 0.0]
+    chord = [2.2, 1.8]
+    theta = [2.0*pi/180, 2.0*pi/180]
+    phi = [0.0, 0.0]
+    ns = 12
+    nc = 1
+    spacing_s = Uniform()
+    spacing_c = Uniform()
+
+    Sref = 30.0
+    cref = 2.0
+    bref = 15.0
+    rref = [0.50, 0.0, 0.0]
+    Vinf = 1.0
+    ref = Reference(Sref, cref, bref, rref, Vinf)
+
+    alpha = 1.0*pi/180
+    beta = 0.0
+    Omega = [0.0; 0.0; 0.0]
+    fs = Freestream(Vinf, alpha, beta, Omega)
+
+    # adjust chord length so x-chord length matches AVL
+    chord = @. chord/cos(theta)
+
+    # vortex rings with symmetry
+    mirror = false
+    symmetric = true
+
+    grid, ratio = wing_to_grid(xle, yle, zle, chord, theta, phi, ns, nc;
+        mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c)
+
+    grids = [grid]
+    ratios = [ratio]
+
+    system = System(grids; ratios)
+
+
+    mydir = @__DIR__
+    savepath = joinpath(mydir, "test_system.jld2")
+    VortexLattice.save_system_to_bson(system, savepath)
+    loaded_system = VortexLattice.load_system_from_bson(savepath)
+    rm(savepath; force=true)
+
+    for name in fieldnames(typeof(system))
+        if name != :sections
+            @test getfield(system, name) == getfield(loaded_system, name)
+        end
+    end
+end
\ No newline at end of file