diff --git a/src/BLUEs.jl b/src/BLUEs.jl
index 70404f3..2d18b79 100644
--- a/src/BLUEs.jl
+++ b/src/BLUEs.jl
@@ -315,16 +315,22 @@ end
         Solve underdetermined problem
 """
 function solve(up::UnderdeterminedProblem)
-
+    if up.y isa DimArray
+            y = UnitfulMatrix(ustrip.(vec(up.y)), unit.(vec(up.y)))
+        else
+            y = up.y
+    end
+   
     if ismissing(up.x₀)
-        n = up.y
+        n = y
     else
         if up.x₀ isa DimArray
-            x₀ = UnitfulMatrix(up.x₀[:])
+            x₀ = UnitfulMatrix(ustrip.(vec(up.x₀)), unit.(vec(up.x₀)))
         else
             x₀ = up.x₀
         end
-        n = up.y - up.E*x₀
+
+        n = y - up.E*x₀
     end
     Cxy = up.Cxx*transpose(up.E)
     Cyy = up.E*Cxy + up.Cnn
@@ -374,7 +380,12 @@ end
     Cost function contribution from observations
 """
 function datacost( x̃::Union{Estimate,DimEstimate}, p::Union{OverdeterminedProblem,UnderdeterminedProblem})
-    n = p.y - p.E*x̃.v
+    if p.y isa DimArray
+        y = UnitfulMatrix(ustrip.(vec(p.y)), unit.(vec(p.y)))
+    else
+        y = p.y
+    end
+    n = y - p.E*x̃.v
     if typeof(p) == UnderdeterminedProblem
         Cnn⁻¹ = inv(p.Cnn) # not possible for NamedTuple
     elseif typeof(p) == OverdeterminedProblem
@@ -388,7 +399,7 @@ end
 """
 function controlcost( x̃::Union{Estimate,DimEstimate}, p::Union{OverdeterminedProblem,UnderdeterminedProblem})
     if p.x₀ isa DimArray
-        Δx = x̃.v - UnitfulMatrix(vec(p.x₀))
+        Δx = x̃.v - UnitfulMatrix(ustrip.(vec(p.x₀)), unit.(vec(p.x₀)))
     else
         Δx = x̃.v - p.x₀
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index f6be1c0..a1d5bec 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -87,7 +87,7 @@ include("test_functions.jl")
     end
 
     @testset "left-uniform problem with prior info" begin
-	y = UnitfulMatrix([-1.9permil])
+	y = UnitfulMatrix([-1.], [permil])
 	σₙ = 0.2permil
 	a = -0.24permil*K^-1
 	γδ = 1.0permil^-2 # to keep units correct, have two γ (tapering) variables
@@ -320,29 +320,31 @@ include("test_functions.jl")
                 @test first(E*UnitfulMatrix(vec(x₀))) .== first(predict(x₀))
             
                 # Julia doesn't have method to make scalar into vector
-                !(y isa AbstractVector) && (y = [y])
+                !(y isa AbstractVector) && (y = UnitfulMatrix(ustrip.([y]), unit.([y])))
 
                 # Does E matrix work properly?
-                ỹ = E*UnitfulMatrix(vec(x))
-                for jj in eachindex(y)
+                xvec = UnitfulMatrix(ustrip.(vec(x)), unit.(vec(x)))
+                ỹ = E*xvec
+                for jj in eachindex(vec(y))
                     @test isapprox(vec(y)[jj],vec(ỹ)[jj])
                 end
-                x̂ = E\y
-                for jj in eachindex(y)
-                    @test isapprox(vec(y)[jj],vec(E*x̂)[jj])
+                yvec = UnitfulMatrix(ustrip.(vec(y)), unit.(vec(y))) #convert from DA
+                x̂ = E\yvec
+                for jj in eachindex(vec(y))
+                    @test isapprox(vec(y)[jj],getindexqty(E*x̂, jj))
                 end
 
                 # now in a position to use BLUEs to solve
                 σₙ = 0.01
                 σₓ = 100.0
 
-                Cnn = Diagonal(fill(σₙ^2,length(y)),unit.(y),unit.(y).^-1)
+                Cnn = Diagonal(fill(σₙ^2,length(y)),unitrange(yvec),unitrange(yvec).^-1)
                 Cxx = Diagonal(fill(σₓ^2,length(x₀)),vec(unit.(x₀)),vec(unit.(x₀)).^-1)
-                problem = UnderdeterminedProblem(UnitfulMatrix(y),E,Cnn,Cxx,x₀)
+                problem = UnderdeterminedProblem(y,E,Cnn,Cxx,x₀)
 
                 # when x₀ is a DimArray, then x̃ is a DimEstimate
                 x̃ = solve(problem) # ::DimEstimate
-                @test within(y[1],vec((E*x̃).v)[1],3σₙ) # within 3-sigma
+                @test within(getindexqty(yvec, 1),vec((E*x̃).v)[1],3σₙ) # within 3-sigma
                 @test cost(x̃,problem) < 5e-2 # no noise in ob
                 @test cost(x̃, problem) == datacost(x̃, problem) + controlcost(x̃, problem)
 
@@ -385,10 +387,12 @@ include("test_functions.jl")
                 y = predict(x)
 
                 # test compatibility
-                @test first(E*UnitfulMatrix(vec(x₀))) .== first(predict(x₀))
+                x₀vec = UnitfulMatrix(ustrip.(vec(x₀)), unit.(vec(x₀)))
+                @test first(E*x₀vec) .== first(predict(x₀))
             
                 # Does E matrix work properly?
-                ỹ = E*UnitfulMatrix(vec(x))
+                xvec = UnitfulMatrix(ustrip.(vec(x)), unit.(vec(x)))
+                ỹ = E*UnitfulMatrix(vec(xvec))
                 for jj in eachindex(y)
                     @test isapprox(vec(y)[jj],vec(ỹ)[jj])
                 end
@@ -405,7 +409,9 @@ include("test_functions.jl")
                 σₓ = 10.0
                 Cnn = Diagonal(fill(σₙ^2,length(y)),vec(unit.(y)),vec(unit.(y).^-1))
                 Cxx = Diagonal(fill(σₓ^2,length(x₀)),vec(unit.(x₀)),vec(unit.(x₀)).^-1)
-                problem = UnderdeterminedProblem(UnitfulMatrix(vec(y)),E,Cnn,Cxx,x₀)
+                
+                yvec = UnitfulMatrix(ustrip.(vec(y)), unit.(vec(y)))
+                problem = UnderdeterminedProblem(yvec,E,Cnn,Cxx,x₀)
 
                 # when x₀ is a DimArray, then x̃ is a DimEstimate
                 x̃ = solve(problem)