diff --git a/pynumdiff/linear_model/_linear_model.py b/pynumdiff/linear_model/_linear_model.py
index 1fd73a9..60cc077 100644
--- a/pynumdiff/linear_model/_linear_model.py
+++ b/pynumdiff/linear_model/_linear_model.py
@@ -6,10 +6,15 @@
 from pynumdiff.finite_difference import first_order as finite_difference
 from pynumdiff.utils import utility
 
+try: import cvxpy
+except ImportError: pass
+
 try:
-    import cvxpy
+    import mosek
+    solver = 'MOSEK' # https://www.mosek.com/
 except ImportError:
-    pass
+    warn("MOSEK not installed, falling back to CVXPY's defaults")
+    solver = None # passing this to solve() allows CVXPY to use whatever it deems best
 
 KERNELS = {'friedrichs': utility._friedrichs_kernel,
            'gaussian': utility._gaussian_kernel}
@@ -299,7 +304,7 @@ def polydiff(x, dt, params=None, options=None, polynomial_order=None, window_siz
 
 
 def __solve_for_A_and_C_given_X_and_Xdot__(X, Xdot, num_integrations, dt, gammaC=1e-1, gammaA=1e-6,
-                                           solver='MOSEK', A_known=None, epsilon=1e-6, rows_of_interest='all'):
+                                           solver=solver, A_known=None, epsilon=1e-6, rows_of_interest='all'):
     """Given state and the derivative, find the system evolution and measurement matrices.
     """
 
@@ -355,7 +360,7 @@ def __integrate_dxdt_hat_matrix__(dxdt_hat, dt):
     return x
 
 
-def _lineardiff(x, dt, N, gamma, solver='MOSEK', weights=None):
+def _lineardiff(x, dt, N, gamma, solver=solver, weights=None):
     """Estimate the parameters for a system xdot = Ax, and use that to calculate the derivative
 
     :param np.array[float] x: time series to differentiate
@@ -405,7 +410,7 @@ def _lineardiff(x, dt, N, gamma, solver='MOSEK', weights=None):
 
 
 def lineardiff(x, dt, params=None, options=None, order=None, gamma=None, window_size=None,
-    sliding=True, step_size=10, kernel='friedrichs', solver='MOSEK'):
+    sliding=True, step_size=10, kernel='friedrichs', solver=solver):
     """Slide a smoothing derivative function across a time series with specified window size.
 
     :param np.array[float] x: array of time series to differentiate
diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py
index def67bc..80a2d10 100644
--- a/pynumdiff/tests/test_diff_methods.py
+++ b/pynumdiff/tests/test_diff_methods.py
@@ -33,7 +33,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
     (first_order, {}), # empty dictionary for the case of no parameters. no params -> no diff in new vs old
     (iterated_first_order, {'num_iterations':5}), (iterated_first_order, [5], {'iterate':True}),
     (second_order, {}),
-    #(lineardiff, {'order':3, 'gamma':5, 'window_size':10, 'solver':'CVXOPT'}),
+    (lineardiff, {'order':3, 'gamma':5, 'window_size':10}), (lineardiff, [3, 5, 10]),
     (polydiff, {'polynomial_order':2, 'window_size':3}), (polydiff, [2, 3]),
     (savgoldiff, {'polynomial_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]),
     (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1]),
@@ -67,7 +67,12 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
                    [(-25, -25), (0, -1), (0, 0), (1, 1)],
                    [(-25, -25), (1, 1), (0, 0), (1, 1)],
                    [(-25, -25), (3, 3), (0, 0), (3, 3)]],
-    #lineardiff: [TBD when #91 is solved],
+    lineardiff: [[(-6, -6), (-5, -6), (-1, -1), (0, 0)],
+                 [(0, 0), (1, 1), (0, 0), (1, 1)],
+                 [(1, 0), (2, 1), (1, 0), (2, 1)],
+                 [(1, 0), (1, 1), (1, 0), (1, 1)],
+                 [(1, 1), (2, 2), (1, 1), (2, 2)],
+                 [(1, 1), (3, 3), (1, 1), (3, 3)]],
     polydiff: [[(-14, -15), (-14, -14), (0, -1), (1, 1)],
                [(-14, -14), (-13, -13), (0, -1), (1, 1)],
                [(-14, -14), (-13, -13), (0, -1), (1, 1)],