diff --git a/docs/src/example/generate-READMEs.jl b/docs/src/example/generate-READMEs.jl
index 073ae7a4..87f4d2c3 100644
--- a/docs/src/example/generate-READMEs.jl
+++ b/docs/src/example/generate-READMEs.jl
@@ -1,6 +1,6 @@
 # find julia files to be converted into README.md
 
-function gernerate_README(examplestoconvert, folder)
+function generate_README(examplestoconvert, folder)
 
     ps = map(examplestoconvert) do str
         fs = filter(s -> split(s, ".")[end] == "jl", readdir(folder*str))
diff --git a/docs/src/example/particles_in_circle/README.md b/docs/src/example/particles_in_circle/README.md
index 1b0ebda4..e44f6ef6 100644
--- a/docs/src/example/particles_in_circle/README.md
+++ b/docs/src/example/particles_in_circle/README.md
@@ -37,7 +37,7 @@ big_particle_simulation = FrequencySimulation([big_particle], source)
 using Plots
 height = 300
 #gr(size=(1.4*height,height))
-pyplot(leg=false, size=(1.4*height,height))
+#pyplot(leg=false, size=(1.4*height,height))
 
 ω = 0.5
 #plot(big_particle_simulation, ω; res=15, bounds = box);
@@ -64,8 +64,8 @@ big_result = run(big_particle_simulation, x, ωs)
 
 #plot(result, lab = "scattering from particles")
 #plot!(big_result,
-#    lab = "scattering from big particle",
-#    title="Compare scattered wave from one big particle, \n and a circle filled with small particles")
+#lab = "scattering from big particle",
+#title="Compare scattered wave from one big particle, \n and a circle filled with small particles")
 ```
 ![The response comparison](plot_response_compare.png)
 
diff --git a/src/physics/acoustics/acoustics.jl b/src/physics/acoustics/acoustics.jl
index b30ec95f..1ad41f13 100644
--- a/src/physics/acoustics/acoustics.jl
+++ b/src/physics/acoustics/acoustics.jl
@@ -146,9 +146,9 @@ function internal_field(x::AbstractVector{T}, p::Particle{Dim,Acoustic{T,Dim}},
         vos = regular_radial_basis(p.medium, ω, order, r)
         us = outgoing_radial_basis(source.medium, ω, order, r)
 
-        internal_coefs = (vs .* (inv(t_mat) * fs) + us .* fs) ./ vos
+        internal_coefs = (vs[:] .* (inv(t_mat) * fs) + us[:] .* fs) ./ vos[:]
         inner_basis = regular_basis_function(p, ω)
 
-        return sum(inner_basis(order, x-origin(p)) .* internal_coefs)
+        return inner_basis(order, x-origin(p)) * internal_coefs
     end
 end
diff --git a/src/physics/acoustics/concentric_capsule.jl b/src/physics/acoustics/concentric_capsule.jl
index 623ded10..0bc92fe7 100644
--- a/src/physics/acoustics/concentric_capsule.jl
+++ b/src/physics/acoustics/concentric_capsule.jl
@@ -66,7 +66,7 @@ function internal_field(x::AbstractArray{T}, p::CapsuleParticle{2,Acoustic{T,2},
         coefs = - q0 * forces .* Ydn.(-Nh:Nh, a0*k1, a0*k1)
         basis = regular_basis_function(p.inner, ω)
 
-        return sum(basis(Nh,x-origin(p)) .* coefs)
+        return basis(Nh,x-origin(p)) * coefs
     else # calculate field in the mid region
         J_coefs = forces .* [
             q0*besselj(n, a0*k0)*diffhankelh1(n, a0*k1) - hankelh1(n, a0*k1)*diffbesselj(n, a0*k0)
@@ -78,8 +78,6 @@ function internal_field(x::AbstractArray{T}, p::CapsuleParticle{2,Acoustic{T,2},
 
         J_basis = regular_basis_function(p.outer, ω)
         H_basis = outgoing_basis_function(p.outer.medium, ω)
-        return sum(
-            J_basis(Nh,x-origin(p)) .* J_coefs .+ H_basis(Nh,x-origin(p)) .* H_coefs
-        )
+        return J_basis(Nh,x-origin(p)) * J_coefs + H_basis(Nh,x-origin(p)) * H_coefs
     end
 end
diff --git a/src/physics/physical_medium.jl b/src/physics/physical_medium.jl
index f97e3029..f9167ea3 100644
--- a/src/physics/physical_medium.jl
+++ b/src/physics/physical_medium.jl
@@ -64,21 +64,21 @@ basislength_to_basisorder(::Type{P},len::Int) where P <: PhysicalMedium{2,1} = I
 
 function outgoing_radial_basis(medium::PhysicalMedium{2,1}, ω::T, order::Integer, r::T) where {T<:Number}
     k = ω/medium.c
-    return hankelh1.(-order:order,k*r)
+    return transpose(hankelh1.(-order:order,k*r))
 end
 
 function outgoing_basis_function(medium::PhysicalMedium{2,1}, ω::T) where {T<:Number}
     return function (order::Integer, x::AbstractVector{T})
         r, θ  = cartesian_to_radial_coordinates(x)
         k = ω/medium.c
-        [hankelh1(m,k*r)*exp(im*θ*m) for m = -order:order]
+        [hankelh1(m,k*r)*exp(im*θ*m) for m = -order:order] |> transpose
     end
 end
 
 function outgoing_radial_basis(medium::PhysicalMedium{3,1}, ω::T, order::Integer, r::T) where {T<:Number}
     k = ω/medium.c
     hs = shankelh1.(0:order,k*r)
-    return  [hs[l+1] for l = 0:order for m = -l:l]
+    return  [hs[l+1] for l = 0:order for m = -l:l] |> transpose
 end
 
 function outgoing_basis_function(medium::PhysicalMedium{3,1}, ω::T) where {T<:Number}
@@ -91,7 +91,7 @@ function outgoing_basis_function(medium::PhysicalMedium{3,1}, ω::T) where {T<:N
 
         lm_to_n = lm_to_spherical_harmonic_index
 
-        return [hs[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l]
+        return [hs[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l] |> transpose
     end
 end
 
@@ -142,7 +142,7 @@ regular_basis_function(medium::PhysicalMedium{Dim,1}, ω::Union{T,Complex{T}}) w
 
 function regular_radial_basis(medium::PhysicalMedium{2,1}, ω::T, order::Integer, r::T) where {T<:Number}
     k = ω / medium.c
-    return besselj.(-order:order,k*r)
+    return transpose(besselj.(-order:order,k*r))
 end
 
 function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::PhysicalMedium{2,1}) where T
@@ -150,7 +150,7 @@ function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::PhysicalMediu
         r, θ  = cartesian_to_radial_coordinates(x)
         k = wavenumber
 
-        return [besselj(m,k*r)*exp(im*θ*m) for m = -order:order]
+        return [besselj(m,k*r)*exp(im*θ*m) for m = -order:order] |> transpose
     end
 end
 
@@ -158,7 +158,7 @@ function regular_radial_basis(medium::PhysicalMedium{3,1}, ω::T, order::Integer
     k = ω / medium.c
     js = sbesselj.(0:order,k*r)
 
-    return [js[l+1] for l = 0:order for m = -l:l]
+    return [js[l+1] for l = 0:order for m = -l:l] |> transpose
 end
 
 function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::PhysicalMedium{3,1}) where T
@@ -170,7 +170,7 @@ function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::PhysicalMediu
 
         lm_to_n = lm_to_spherical_harmonic_index
 
-        return [js[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l]
+        return [js[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l] |> transpose
     end
 end
 
diff --git a/src/physics/special_functions.jl b/src/physics/special_functions.jl
index e641290a..435bcbfa 100644
--- a/src/physics/special_functions.jl
+++ b/src/physics/special_functions.jl
@@ -54,30 +54,30 @@ function diffsbessel(f::Function,n::Number,z::Number)
     return f(n-1,z) - (n+1) * f(n,z) / z
 end
 
-"""
-    diffsbessel(f::Function,m,x,n::Int)
+# """
+#     diffsbessel(f::Function,m,x,n::Int)
 
-Differentiates 'n' times any spherical Bessel function 'f' of order 'm' and at the argument 'x'. Uses the formula 
+# Differentiates 'n' times any spherical Bessel function 'f' of order 'm' and at the argument 'x'. Uses the formula 
     
-``(1/z d/dz)^n (z^{m+1} f_m(z)) = z^{m-n+1} f_{m-n}``
+# ``(1/z d/dz)^n (z^{m+1} f_m(z)) = z^{m-n+1} f_{m-n}``
 
-which leads to     
+# which leads to     
 
-``(1/z d/dz)^{n} (z^{m+1} f_m(z)) = (1/z d/dz)^{n-1} (1/z d/dz) (z^{m+1} f_m(z))``
+# ``(1/z d/dz)^{n} (z^{m+1} f_m(z)) = (1/z d/dz)^{n-1} (1/z d/dz) (z^{m+1} f_m(z))``
     
-"""
-function diffsbessel(f::Function,m::Number,z,n::Int)
-    if n == 0
-        return f(m, z)
-    elseif n > 0
-        n = n - 1
-        return diffsbessel(f,m,z)
+# """
+# function diffsbessel(f::Function,m::Number,z,n::Int)
+#     if n == 0
+#         return f(m, z)
+#     elseif n > 0
+#         n = n - 1
+#         return diffsbessel(f,m,z)
 
         
-    else
-        error("Can not differentiate a negative number of times")
-    end
-end
+#     else
+#         error("Can not differentiate a negative number of times")
+#     end
+# end
 
 function diffsbesselj(n::Number,z::Number)
     return if n == 0
diff --git a/src/simulation.jl b/src/simulation.jl
index 4f919f67..6e85a2c5 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -244,7 +244,7 @@ function field(sim::FrequencySimulation{Dim,P}, ω::Number, x_vec::Vector{V}, a_
     function sum_basis(x)
         sum(eachindex(sim.particles)) do i
             p = sim.particles[i]
-            sum(a_vec[:,i] .* basis(N, x-origin(p)))
+            basis(N, x-origin(p)) * a_vec[:,i]
         end
     end
     map(x_vec) do x
diff --git a/src/source.jl b/src/source.jl
index d40d72e6..3492b604 100644
--- a/src/source.jl
+++ b/src/source.jl
@@ -141,7 +141,7 @@ function regular_spherical_source(medium::PhysicalMedium{Dim},regular_coefficien
 
     function source_field(x,ω)
         vs = regular_basis_function(medium, ω)(coeff_order,x-position)
-        return amplitude * sum(vs .* regular_coefficients)
+        return amplitude * (vs * regular_coefficients)
 
     end
 
@@ -171,7 +171,8 @@ function source_expand(source::AbstractSource, centre::AbstractVector{T}; basis_
         vs = regular_basis_function(source.medium, ω)
         regular_coefficients = regular_spherical_coefficients(source)
 
-        res = sum(regular_coefficients(basis_order,centre,ω) .* vs(basis_order, x - centre), dims = 1)[:]
+        res = vs(basis_order, x - centre) * regular_coefficients(basis_order,centre,ω)
+        # res = sum(regular_coefficients(basis_order,centre,ω) .* vs(basis_order, x - centre), dims = 1)[:]
         if length(res) == 1
             res = res[1]
         end
diff --git a/test/acoustic_physics_tests.jl b/test/acoustic_physics_tests.jl
index bae04739..46eedaae 100644
--- a/test/acoustic_physics_tests.jl
+++ b/test/acoustic_physics_tests.jl
@@ -49,13 +49,13 @@ end
     vs = regular_basis_function(medium, ω)(4*order,r)
     us = outgoing_basis_function(medium, ω)(order,r + d)
 
-    @test maximum(abs.(U * vs - us) ./ abs.(us)) < 3e-6 # for linux: 4e-7
+    @test maximum(abs.(U * vs[:] - us[:]) ./ abs.(us[:])) < 3e-6 # for linux: 4e-7
 
     V = regular_translation_matrix(medium, 4*order, order, ω, d)
     v1s = regular_basis_function(medium, ω)(4*order,r)
     v2s = regular_basis_function(medium, ω)(order,r + d)
 
-    @test maximum(abs.(V * v1s - v2s) ./ abs.(v2s)) < 1e-8 # for linux: 4e-7
+    @test maximum(abs.(V * v1s[:] - v2s[:]) ./ abs.(v2s[:])) < 1e-8 # for linux: 4e-7
 
     # Test 2D outgoing translation matrix
     ω = rand() + 0.1
@@ -71,13 +71,13 @@ end
     vs = regular_basis_function(medium, ω)(larger_order,r)
     us = outgoing_basis_function(medium, ω)(order,r + d)
 
-    @test maximum(abs.(U * vs - us) ./ abs.(us)) < 1e-9
+    @test maximum(abs.(U * vs[:] - us[:]) ./ abs.(us[:])) < 1e-9
 
     V = regular_translation_matrix(medium, larger_order, order, ω, d)
     v1s = regular_basis_function(medium, ω)(larger_order,r)
     v2s = regular_basis_function(medium, ω)(order,r + d)
 
-    @test maximum(abs.(V * v1s - v2s) ./ abs.(v2s)) < 1e-10
+    @test maximum(abs.(V * v1s[:] - v2s[:]) ./ abs.(v2s[:])) < 1e-10
 
 end
 
diff --git a/test/run_examples.jl b/test/run_examples.jl
index 8145bfc9..8de94971 100644
--- a/test/run_examples.jl
+++ b/test/run_examples.jl
@@ -13,7 +13,7 @@ using MultipleScattering, Test
 
     include(folder*"generate-READMEs.jl")
 
-    gernerate_README(examplestoconvert, folder)
+    generate_README(examplestoconvert, folder)
 
     @test true
 end