diff --git a/Project.toml b/Project.toml
index 0375150..4a4e901 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,15 +1,37 @@
 name = "DualArrays"
 uuid = "429e4e16-f749-45f3-beec-30742fae38ce"
+version = "0.2.0"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.1.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[compat]
+ChainRules = "1.72.6"
+ChainRulesCore = "1.26.0"
+DifferentialEquations = "7.17.0"
+ForwardDiff = "1.2.2"
+Plots = "1.41.1"
+SparseArrays = "1.10"
+
+[extras]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[targets]
+dev = ["Revise"]
+examples = ["Plots", "DifferentialEquations", "ForwardDiff", "BandedMatrices"]
+test = ["Test", "ForwardDiff", "BandedMatrices", "SparseArrays"]
diff --git a/README.md b/README.md
index 8cc9eeb..078a9f9 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,19 @@
 # DualArrays.jl
- A package for working with arrays of dual numbers with sparsity
+DualArrays.jl is a package that provides the `DualVector` structure for use in forward mode automatic differentiation (autodiff). Existing forward mode autodiff implementations such as `ForwardDiff.jl` make use of a `Dual` structure with a vector of dual parts.
+
+There are some limitations of this when differentiating vector valued functions, as the dual components of each element will be each treated as a separate (dense) vector rather than a Jacobian. This misses the opportunity to exploit sparse Jacobian structures such as banded structures that appear when solving systems of ODEs. 
+
+`DualArrays.jl` provides a new structure, `DualVector`, consisting of a vector of real parts and a jacobian that can be any matrix structure in the Julia ecosystem. This carries over many of the optimisations provided by sparse matrix structures to the forward-mode autodiff process. The package also comes with its own `Dual` type for elementwise indexing or vector -> scalar functions. An efficient implementation of differentiating a function might be as follows:
+
+```julia
+using FillArrays
+
+function gradient(f::Function, x::Vector)
+    dx = DualVector(x, Eye(length(x)))
+    return f(dx).jacobian
+end
+```
+
+See the examples folder for more use cases.
+
+Differentiation rules are mostly provided by the `ChainRules.jl` autodiff backend.
diff --git a/examples/newtonpendulum.jl b/examples/newtonpendulum.jl
index 81d727b..7802ef8 100644
--- a/examples/newtonpendulum.jl
+++ b/examples/newtonpendulum.jl
@@ -48,8 +48,8 @@ end
 x0 = zeros(Float64, N)
 
 #Solve and plot both solution and LHS ('deviation' from system)
-@time sol = newton_method_forwarddiff(f, x0, 100)
-@time sol = newton_method_dualvector(f, x0, 100)
+@time sol = newton_method_forwarddiff(f, x0, 10)
+@time sol = newton_method_dualvector(f, x0, 10)
 plot(0:ts:Tmax, [a; sol; b])
 plot!(0:ts:Tmax, [0; f(sol); 0])
 
diff --git a/src/DualArrays.jl b/src/DualArrays.jl
index e4074d5..7feb3d5 100644
--- a/src/DualArrays.jl
+++ b/src/DualArrays.jl
@@ -1,103 +1,25 @@
-module DualArrays
-export DualVector
-
-import Base: +, ==, getindex, size, broadcast, axes, broadcasted, show, sum,
-             vcat, convert, *
-using LinearAlgebra, ArrayLayouts, BandedMatrices, FillArrays
-
-struct Dual{T, Partials <: AbstractVector{T}} <: Real
-    value::T
-    partials::Partials
-end
-
-==(a::Dual, b::Dual) = a.value == b.value && a.partials == b.partials
-
-sparse_getindex(a...) = layout_getindex(a...)
-sparse_getindex(D::Diagonal, k::Integer, ::Colon) = OneElement(D.diag[k], k, size(D,2))
-sparse_getindex(D::Diagonal, ::Colon, j::Integer) = OneElement(D.diag[j], j, size(D,1))
-
 """
-reprents a vector of duals given by
-    
-    values + jacobian * [ε_1,…,ε_n].
-
-For now the entries just return the values.
-"""
-
-struct DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
-    value::Vector{T}
-    jacobian::M
-    function DualVector(value::Vector{T},jacobian::M) where {T, M <: AbstractMatrix{T}}
-        if(size(jacobian)[1] != length(value))
-            x,y = length(value),size(jacobian)[1]
-            throw(ArgumentError("vector length must match number of rows in jacobian.\n
-            vector length: $x \n
-            no. of jacobian rows: $y"))
-        end
-        new{T,M}(value,jacobian)
-    end
-end
-
-function DualVector(value::AbstractVector, jacobian::AbstractMatrix)
-    T = promote_type(eltype(value), eltype(jacobian))
-    DualVector(convert(Vector{T}, value), convert(AbstractMatrix{T}, jacobian))
-end
+    DualArrays
 
+A Julia package for efficient automatic differentiation using dual numbers and dual arrays.
 
+This package provides:
+- `Dual`: A dual number type for storing values and their derivatives
+- `DualVector`: A vector of dual numbers represented with a Jacobian matrix
 
-function getindex(x::DualVector, y::Int)
-    Dual(x.value[y], sparse_getindex(x.jacobian,y,:))
-end
-
-function getindex(x::DualVector, y::UnitRange)
-    newval = x.value[y]
-    newjac = sparse_getindex(x.jacobian,y,:)
-    DualVector(newval, newjac)
-end
-size(x::DualVector) = length(x.value)
-axes(x::DualVector) = axes(x.value)
-+(x::DualVector,y::DualVector) = DualVector(x.value + y.value, x.jacobian + y.jacobian)
-*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
-
-broadcasted(::typeof(sin),x::DualVector) = DualVector(sin.(x.value),Diagonal(cos.(x.value))*x.jacobian)
-
-function broadcasted(::typeof(*),x::DualVector,y::DualVector)
-    newval = x.value .* y.value
-    newjac = x.value .* y.jacobian + y.value .* x.jacobian
-    DualVector(newval,newjac)
-end
-
-function sum(x::DualVector)
-    n = length(x.value)
-    Dual(sum(x.value), vec(sum(x.jacobian; dims=1)))
-end
-
-_jacobian(d::Dual) = permutedims(d.partials)
-_jacobian(d::DualVector, ::Int) = d.jacobian
-_jacobian(x::Number, N::Int) = zeros(typeof(x), 1, N)
+Differentiation rules are mostly provided by ChainRules.jl.
+"""
+module DualArrays
 
-_value(d::DualVector) = d.value
-_value(x::Number) = x
+export DualVector, Dual
 
-function vcat(x::Union{Dual, DualVector}...)
-    if length(x) == 1
-        return x[1]
-    end
-    value = vcat((d.value for d in x)...)
-    jacobian = vcat((_jacobian(d) for d in x)...)
-    DualVector(value,jacobian)
-end
+import Base: +, ==, getindex, size, axes, broadcasted, show, sum, vcat, convert, *
 
-_size(x::Real) = 1
-_size(x::DualVector) = size(x)
+using LinearAlgebra, ArrayLayouts, FillArrays
 
-function vcat(x::Union{Real, DualVector}...)
-    cols = max((_size(i) for i in x)...)
-    val = vcat((_value(i) for i in x)...)
-    jac = vcat((_jacobian(i, cols) for i in x)...)
-    DualVector(val, jac)
+include("types.jl")
+include("indexing.jl")
+include("arithmetic.jl")
+include("utilities.jl")
 end
 
-show(io::IO,::MIME"text/plain", x::DualVector) = (print(io,x.value); print(io," + "); print(io,x.jacobian);print("𝛜"))
-end
-# module DualArrays
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
new file mode 100644
index 0000000..1fb94f7
--- /dev/null
+++ b/src/arithmetic.jl
@@ -0,0 +1,170 @@
+# Arithmetic operations for DualArrays.jl
+
+using LinearAlgebra, ChainRules
+
+"""
+Addition of DualVectors.
+"""
+Base.:+(x::DualVector, y::DualVector) = DualVector(x.value + y.value, x.jacobian + y.jacobian)
+Base.:+(x::DualVector, y::AbstractVector) = DualVector(x.value + y, x.jacobian)
+Base.:+(x::AbstractVector, y::DualVector) = DualVector(x + y.value, y.jacobian)
+
+"""
+Matrix multiplication with a DualVector.
+"""
+Base.:*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
+
+"""
+Broadcasted multiplication of two DualVectors using the product rule.
+"""
+function Base.broadcasted(::typeof(*), x::DualVector, y::DualVector)
+    newval = x.value .* y.value
+    newjac = Diagonal(x.value) * y.jacobian + Diagonal(y.value) * x.jacobian
+    DualVector(newval, newjac)
+end
+
+"""
+
+this section attempts to define broadcasting rules on DualVectors for functions
+that either:
+
+- take a single real argument (function applied to each element)
+- are binary operations (binary operation applied to scalar and each element)
+
+This implementation is loosely based on
+https://juliadiff.org/ChainRulesOverloadGeneration.jl/dev/examples/forward_mode.html
+
+"""
+
+"""
+Defines how a frule in ChainRules.jl for a scalar function f should be broadcasted over a DualVector.
+"""
+function broadcast_rule(f, d::DualVector)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+
+    @inbounds for (i, x) in pairs(d.value)
+        y, dy = ChainRules.frule((ChainRules.NoTangent(), d.jacobian[i, :]), f, x)
+        val[i] = y
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+"""
+Defines how a frule in ChainRules.jl for a binary
+operation a (+) b on reals a, b should be broadcasted over:
+- a (+) d (d a DualVector)
+- d (+) a (d a DualVector)
+"""
+
+function broadcast_rule(f, d::DualVector, x::Real)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+    z = zero(jac[1, :])
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), d.jacobian[i, :], z),
+             f, y, x)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+function broadcast_rule(f, x::Real, d::DualVector)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+    z = zero(jac[1, :])
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), z, d.jacobian[i, :]),
+             f, x, y)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+"""
+Extend for broadcasting Dual and DualVector
+"""
+
+function broadcast_rule(f, d::DualVector, x::Dual)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+    z = x.partials
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), d.jacobian[i, :], z),
+             f, y, x.value)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+function broadcast_rule(f, x::Dual, d::DualVector)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+    z = x.partials
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), z, d.jacobian[i, :]),
+             f, x.value, y)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+# a set of defined broadcasts to avoid duplicate definitions
+defined = Set{DataType}()
+
+# Get all applicable frules defined in ChainRules
+# and define broadcasted versions for DualVector using vector_rule
+
+frules = methods(ChainRules.frule)
+for f in frules
+    # get signatures for each function with a frule
+    sig = Base.unwrap_unionall(Base.tuple_type_tail(f.sig))
+    
+    # split into operation and args, filter for
+    # single argument functions and binary operations that can act on real numbers
+    op, args = sig.parameters[2], sig.parameters[3:end]
+
+    isconcretetype(op) || continue
+    op in defined && continue
+
+    # if it is a single argument function...
+    if length(args) == 1
+        args[1] isa Type || continue
+        Real <: args[1] || continue
+
+        @eval Base.broadcasted(fn::$op, d::DualVector) = broadcast_rule(fn, d)
+        push!(defined, op)
+    end
+
+    # if it is a binary operation...
+    if length(args) == 2
+        args[1] isa Type || continue
+        args[2] isa Type || continue
+        Real <: args[1] && Real <: args[2] || continue
+
+        
+        @eval Base.broadcasted(fn::$op, d::DualVector, y::Real) = broadcast_rule(fn, d, y)
+        @eval Base.broadcasted(fn::$op, y::Real, d::DualVector) = broadcast_rule(fn, y, d)
+        @eval Base.broadcasted(fn::$op, d::DualVector, du::Dual) = broadcast_rule(fn, d, du)
+        @eval Base.broadcasted(fn::$op, du::Dual, d::DualVector) = broadcast_rule(fn, du, d)
+        push!(defined, op)
+    end
+end
\ No newline at end of file
diff --git a/src/indexing.jl b/src/indexing.jl
new file mode 100644
index 0000000..7036091
--- /dev/null
+++ b/src/indexing.jl
@@ -0,0 +1,33 @@
+# Indexing operations for DualArrays.jl
+
+using ArrayLayouts, FillArrays
+
+sparse_getindex(a...) = layout_getindex(a...)
+sparse_getindex(D::Diagonal, k::Integer, ::Colon) = OneElement(D.diag[k], k, size(D, 2))
+sparse_getindex(D::Diagonal, ::Colon, j::Integer) = OneElement(D.diag[j], j, size(D, 1))
+
+"""
+Extract a single Dual number from a DualVector at position y.
+"""
+function Base.getindex(x::DualVector, y::Int)
+    Dual(x.value[y], sparse_getindex(x.jacobian, y, :))
+end
+
+"""
+Extract a sub-DualVector from a DualVector using a range.
+"""
+function Base.getindex(x::DualVector, y::UnitRange)
+    newval = x.value[y]
+    newjac = sparse_getindex(x.jacobian, y, :)
+    DualVector(newval, newjac)
+end
+
+"""
+Return the size of the DualVector (length of the value vector).
+"""
+Base.size(x::DualVector) = (length(x.value),)
+
+"""
+Return the axes of the DualVector.
+"""
+Base.axes(x::DualVector) = axes(x.value)
\ No newline at end of file
diff --git a/src/types.jl b/src/types.jl
new file mode 100644
index 0000000..144d713
--- /dev/null
+++ b/src/types.jl
@@ -0,0 +1,59 @@
+# Core type definitions for DualArrays.jl
+
+"""
+    Dual{T, Partials <: AbstractVector{T}} <: Real
+
+A dual number type that stores a value and its partials (derivatives).
+
+# Fields
+- `value::T`: The primal value
+- `partials::Partials`: The partial derivatives as a vector
+"""
+struct Dual{T, Partials <: AbstractVector{T}} <: Real
+    value::T
+    partials::Partials
+end
+
+"""
+    DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
+
+Represents a vector of dual numbers given by:
+    
+    values + jacobian * [ε₁, …, εₙ]
+
+For now the entries just return the values when indexed.
+
+# Fields
+- `value::Vector{T}`: The primal values
+- `jacobian::M`: The Jacobian matrix containing partial derivatives
+
+# Constructor
+    DualVector(value::Vector{T}, jacobian::M) where {T, M <: AbstractMatrix{T}}
+
+Constructs a DualVector, ensuring that the vector length matches the number of rows in the Jacobian.
+"""
+struct DualVector{T, V <: AbstractVector{T},M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
+    value::V
+    jacobian::M
+
+    function DualVector(value::V, jacobian::M) where {T, V <: AbstractVector{T}, M <: AbstractMatrix{T}}
+        if size(jacobian, 1) != length(value)
+            x, y = length(value), size(jacobian, 1)
+            throw(ArgumentError("vector length must match number of rows in jacobian.\n" *
+                               "vector length: $x\n" *
+                               "no. of jacobian rows: $y"))
+        end
+        new{T,V, M}(value, jacobian)
+    end
+end
+
+"""
+Constructor that forces type compatibility
+"""
+function DualVector(value::AbstractVector, jacobian::AbstractMatrix)
+    T = promote_type(eltype(value), eltype(jacobian))
+    DualVector(convert(Vector{T}, value), convert(AbstractMatrix{T}, jacobian))
+end
+
+# Basic equality for Dual numbers
+==(a::Dual, b::Dual) = a.value == b.value && a.partials == b.partials
\ No newline at end of file
diff --git a/src/utilities.jl b/src/utilities.jl
new file mode 100644
index 0000000..44818c3
--- /dev/null
+++ b/src/utilities.jl
@@ -0,0 +1,56 @@
+# Miscellaneous functions for DualArrays.jl
+
+"""
+Sum all elements of a DualVector, returning a single Dual number.
+"""
+function Base.sum(x::DualVector)
+    n = length(x.value)
+    Dual(sum(x.value), vec(sum(x.jacobian; dims=1)))
+end
+
+# Helper functions for vcat operations
+_jacobian(d::Dual) = permutedims(d.partials)
+_jacobian(d::DualVector) = d.jacobian
+_jacobian(d::DualVector, ::Int) = d.jacobian
+_jacobian(x::Number, N::Int) = zeros(typeof(x), 1, N)
+
+_value(v::AbstractVector) = v
+_value(d::DualVector) = d.value
+_value(d::Dual) = d.value
+_value(x::Number) = x
+
+_size(x::Real) = 1
+_size(x::DualVector) = length(x.value)
+
+"""
+Vertically concatenate Dual numbers and DualVectors.
+"""
+function Base.vcat(x::Union{Dual, DualVector}...)
+    if length(x) == 1
+        return x[1]
+    end
+    value = vcat((d.value for d in x)...)
+    jacobian = vcat((_jacobian(d) for d in x)...)
+    DualVector(value, jacobian)
+end
+
+"""
+Vertically concatenate Real numbers and DualVectors.
+"""
+
+function Base.vcat(x::Union{Real, DualVector}...)
+    # Avoid stack overflow
+    if all(i -> i isa Real, x)
+        return Base.cat(x..., dims=1)
+    end
+    cols = maximum((_size(i) for i in x))
+    val = vcat((_value(i) for i in x)...)
+    jac = vcat((_jacobian(i, cols) for i in x)...)
+    DualVector(val, jac)
+end
+
+"""
+Custom display method for DualVectors.
+"""
+Base.show(io::IO, ::MIME"text/plain", x::DualVector) = 
+    (print(io, x.value); print(io, " + "); print(io, x.jacobian); print(io, "𝛜"))
\ No newline at end of file
diff --git a/test/broadcast_test.jl b/test/broadcast_test.jl
new file mode 100644
index 0000000..df0062b
--- /dev/null
+++ b/test/broadcast_test.jl
@@ -0,0 +1,74 @@
+using DualArrays, Test, SparseArrays, LinearAlgebra
+
+@testset "Broadcasting" begin
+    d = DualVector([1.0, 2.0, 3.0], [1 0 0; 0 1 0; 0 0 1])
+
+    # Test broadcasting with single argument functions
+    s = sin.(d)
+    c = cos.(d)
+    e = exp.(d)
+    l = log.(d)
+    a = tanh.(d)
+
+    @test s isa DualVector
+    @test c isa DualVector
+    @test e isa DualVector
+    @test l isa DualVector
+    @test a isa DualVector
+
+    @test s.value ≈ sin.([1.0, 2.0, 3.0])
+    @test c.value ≈ cos.([1.0, 2.0, 3.0])
+    @test e.value ≈ exp.([1.0, 2.0, 3.0])
+    @test l.value ≈ log.([1.0, 2.0, 3.0])
+    @test a.value ≈ tanh.([1.0, 2.0, 3.0])
+
+    @test s.jacobian ≈ Diagonal(cos.([1.0, 2.0, 3.0]))
+    @test c.jacobian ≈ Diagonal(-sin.([1.0, 2.0, 3.0]))
+    @test e.jacobian ≈ Diagonal(exp.([1.0, 2.0, 3.0]))
+    @test l.jacobian ≈ Diagonal(1.0 ./ [1.0, 2.0, 3.0])
+    @test a.jacobian ≈ Diagonal(1.0 .- tanh.([1.0, 2.0, 3.0]).^2)
+
+    # Test broadcasting with binary operations
+    a = d .+ 2.0
+    s = 3.0 .- d
+    m = d .* 4.0
+    div = 8.0 ./ d
+
+    @test a isa DualVector
+    @test s isa DualVector 
+    @test m isa DualVector
+    @test div isa DualVector
+
+    @test a.value ≈ [3.0, 4.0, 5.0]
+    @test s.value ≈ [2.0, 1.0, 0.0]
+    @test m.value ≈ [4.0, 8.0, 12.0]
+    @test div.value ≈ [8.0, 4.0, 8/3.0]
+
+    @test a.jacobian ≈ d.jacobian
+    @test s.jacobian ≈ -d.jacobian
+    @test m.jacobian ≈ 4.0 * d.jacobian
+    @test div.jacobian ≈ -8.0 ./ (d.value .^ 2) .* d.jacobian
+
+    # Test broadcasting with binary operations on a Dual and DualVector
+    x = Dual(2.0, [1.0, 2.0, 3.0])
+
+    a = x .+ d
+    s = d .- x
+    m = x .* d
+    div = d ./ x
+
+    @test a isa DualVector
+    @test s isa DualVector
+    @test m isa DualVector
+    @test div isa DualVector
+
+    @test a.value ≈ [3.0, 4.0, 5.0]
+    @test s.value ≈ [-1.0, 0.0, 1.0]
+    @test m.value ≈ [2.0, 4.0, 6.0]
+    @test div.value ≈ [0.5, 1.0, 1.5]
+
+    @test a.jacobian ≈ [2.0 2.0 3.0; 1.0 3.0 3.0; 1.0 2.0 4.0]
+    @test s.jacobian ≈ [0.0 -2.0 -3.0; -1.0 -1.0 -3.0; -1.0 -2.0 -2.0]
+    @test m.jacobian ≈ [3.0 2.0 3.0; 2.0 6.0 6.0; 3.0 6.0 11.0]
+    @test div.jacobian ≈ [0.25 -0.5 -0.75; -0.5 -0.5 -1.5; -0.75 -1.5 -1.75]
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index d1cba3c..12a6fc9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,41 +2,53 @@ using DualArrays, Test, LinearAlgebra, ForwardDiff, BandedMatrices
 using DualArrays: Dual
 
 @testset "DualArrays" begin
-    @test_throws ArgumentError DualVector([1,2],I(3))
     
-    v = DualVector([1,2, 3], [1 2 3; 4 5 6;7 8 9])
-    @test v[1] isa Dual
-    @test v[1] == Dual(1,[1,2,3])
-    @test v[2] == Dual(2,[4,5,6])
-    @test v[3] == Dual(3,[7,8,9])
-    @test_throws BoundsError v[4]
-    @test v == DualVector([1,2, 3], [1 2 3; 4 5 6;7 8 9])
+    @testset "Type Definition" begin
+        @test_throws ArgumentError DualVector([1,2],I(3))
+    end
+    
+    @testset "Indexing" begin
+        v = DualVector([1, 2, 3], [1 2 3; 4 5 6;7 8 9])
+        @test v[1] isa Dual
+        @test v[1] == Dual(1,[1,2,3])
+        @test v[2] == Dual(2,[4,5,6])
+        @test v[3] == Dual(3,[7,8,9])
+        @test_throws BoundsError v[4]
+        @test v == DualVector([1,2, 3], [1 2 3; 4 5 6;7 8 9])
 
-    w = v + v
-    @test w == DualVector([2,4,6],[2 4 6;8 10 12;14 16 18])
-    @test w.jacobian == 2v.jacobian
+        x,y = v[1:2],v[2:3]
+        @test x == DualVector([1,2],[1 2 3;4 5 6])
+        @test y == DualVector([2,3],[4 5 6;7 8 9])
 
-    @test sin.(v) isa DualVector
-    @test sin.(v).value == [sin(1), sin(2), sin(3)]
-    @test sin.(v).jacobian == Diagonal(cos.(v.value)) * v.jacobian
+        n = 10
+        v = DualVector(1:n, I(n))
+        @test v[2:end].jacobian isa BandedMatrix
 
-    x,y = v[1:2],v[2:3]
-    @test x == DualVector([1,2],[1 2 3;4 5 6])
-    @test y == DualVector([2,3],[4 5 6;7 8 9])
-    @test x .* y == DualVector([2,6],[6 9 12;26 31 36])
+        @test sum(v[1:end-1] .* v[2:end]).partials == ForwardDiff.gradient(v -> sum(v[1:end-1] .* v[2:end]), 1:n)
+    end
     
-    @test sum(x .* y) isa Dual
-    @test sum(x .* y) == Dual(8,[32,40,48])
-
-    n = 10
-    v = DualVector(1:n, I(n))
-    @test v[2:end].jacobian isa BandedMatrix
+    @testset "Arithmetic" begin
+        v = DualVector([1, 2, 3], [1 2 3; 4 5 6;7 8 9])
+        w = v + v
+        @test w == DualVector([2,4,6],[2 4 6;8 10 12;14 16 18])
+        @test w.jacobian == 2v.jacobian
 
-    @test sum(v[1:end-1] .* v[2:end]).partials == ForwardDiff.gradient(v -> sum(v[1:end-1] .* v[2:end]), 1:n)
+        x = Dual(1, [1, 2, 3])
+        y = DualVector([2, 3], [4 5 6;7 8 9])
 
-    x = Dual(1, [1, 2, 3])
-    y = DualVector([2, 3], [4 5 6;7 8 9])
-    @test vcat(x) == x
-    @test vcat(x, x) == DualVector([1, 1], [1 2 3;1 2 3])
-    @test vcat(x, y) == DualVector([1, 2, 3], [1 2 3;4 5 6;7 8 9])
+        @test x .* y == DualVector([2,3],[6 9 12;10 14 18])
+        
+        @test sum(x .* y) isa Dual
+        @test sum(x .* y) == Dual(5,[16,23,30])
+    end
+    
+    @testset "vcat" begin
+        x = Dual(1, [1, 2, 3])
+        y = DualVector([2, 3], [4 5 6;7 8 9])
+        @test vcat(x) == x
+        @test vcat(x, x) == DualVector([1, 1], [1 2 3;1 2 3])
+        @test vcat(x, y) == DualVector([1, 2, 3], [1 2 3;4 5 6;7 8 9])
+    end
+    
+    include("broadcast_test.jl")
 end
\ No newline at end of file