diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index bade8e8..d36be0c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -15,7 +15,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
           - 'nightly'
         os:
diff --git a/.travis.yml b/.travis.yml
index 14cac2a..53a3d27 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,9 +4,8 @@ os:
   - linux
   - osx
 julia:
-  - 1.0
-  - 1.1
-  - 1.2
+  - 1.6
+  - 1.7
   - nightly
 matrix:
   allow_failures:
diff --git a/Project.toml b/Project.toml
index ea3903b..40d59b0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,16 +7,17 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 
 [compat]
 Combinatorics = "1"
 Distances = "0.9,0.10"
 NearestNeighbors = "0.4"
-julia = "1"
+LinearSolve = "1"
+julia = "1.6"
 
 [extras]
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random"]
+test = ["Test"]
diff --git a/src/ScatteredInterpolation.jl b/src/ScatteredInterpolation.jl
index 86c27fa..bb64138 100755
--- a/src/ScatteredInterpolation.jl
+++ b/src/ScatteredInterpolation.jl
@@ -1,6 +1,6 @@
 module ScatteredInterpolation
 
-using Distances, NearestNeighbors, Combinatorics, LinearAlgebra
+using Distances, NearestNeighbors, Combinatorics, LinearAlgebra, LinearSolve
 
 export interpolate,
     evaluate
diff --git a/src/rbf.jl b/src/rbf.jl
index f956350..3898430 100755
--- a/src/rbf.jl
+++ b/src/rbf.jl
@@ -185,7 +185,7 @@ end
 function interpolate(rbf::Union{T, AbstractVector{T}} where T <: AbstractRadialBasisFunction,
                      points::AbstractArray{<:Real,2},
                      samples::AbstractArray{<:Number,N};
-                     metric = Euclidean(), returnRBFmatrix::Bool = false,
+                     metric = Euclidean(), returnRBFmatrix::Bool = false, linsolve = nothing,
                      smooth::Union{S, AbstractVector{S}} = false) where {N} where {S<:Number}
 
     #hinder smooth from being set to true and interpreted as the value 1 
@@ -198,7 +198,7 @@ function interpolate(rbf::Union{T, AbstractVector{T}} where T <: AbstractRadialB
     A = evaluateRBF!(A, rbf, smooth)
 
     # Solve for the weights
-    itp = solveForWeights(A, points, samples, rbf, metric)
+    itp = solveForWeights(A, points, samples, rbf, metric; linsolve)
 
     # Create and return an interpolation object
     if returnRBFmatrix    # Return matrix A
@@ -245,13 +245,15 @@ end
 
 @inline function solveForWeights(A, points, samples,
                                     rbf::Union{T, AbstractVector{T}} where T <: RadialBasisFunction,
-                                    metric)
-    w = A\samples
-    RBFInterpolant(w, points, rbf, metric)
+                                    metric; linsolve)
+    prob = LinearProblem(A, samples)
+    sol = solve(prob, linsolve)
+
+    RBFInterpolant(sol.u, points, rbf, metric)
 end
 @inline function solveForWeights(A, points, samples,
                                     rbf::Union{T, AbstractVector{T}} where T <: Union{GeneralizedRadialBasisFunction, RadialBasisFunction},
-                                    metric)
+                                    metric; linsolve)
     # Use the maximum degree among the generalized RBF:s
     P = getPolynomial(rbf, points)
 
diff --git a/test/rbf.jl b/test/rbf.jl
index 0814748..8c7fbe3 100755
--- a/test/rbf.jl
+++ b/test/rbf.jl
@@ -39,15 +39,18 @@ radialBasisFunctions = (Gaussian(2),
 
         @test r.(data) ≈ f.(data)
 
-        for points in (arrayPoints, adjointPoints)
-            itp = interpolate(r, points, data)
-
-            # Check that we get back the original data at the sample points and that we get
-            # close when evaluating near the sampling points
-            ev = evaluate(itp, points)
-            @test ev ≈ data
-            ev = evaluate(itp, points .+ 0.001*randn(size(points)))
-            @test all(isapprox.(data, ev, atol = 1e-2))
+        offset = 2e-4
+        for linsolve in (nothing, IterativeSolversJL_GMRES())
+            for points in (arrayPoints, adjointPoints)
+                itp = interpolate(r, points, data; linsolve = linsolve)
+
+                # Check that we get back the original data at the sample points and that we get
+                # close when evaluating near the sampling points
+                ev = evaluate(itp, points)
+                @test ev ≈ data
+                ev = evaluate(itp, points .+ offset)
+                @test all(isapprox.(data, ev, atol = 1e-2))
+            end
         end
     end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 76b70c7..04d311f 100755
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,12 +1,4 @@
-using ScatteredInterpolation, Random
-
-Random.seed!(2)
-
-if VERSION < v"0.7-"
-    using Base.Test
-else
-    using Test
-end
+using ScatteredInterpolation, Test, LinearSolve
 
 include("rbf.jl")
 include("idw.jl")