diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index c34fd24885..fb65274828 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1291,6 +1291,8 @@ function linearization_function(sys::AbstractSystem, inputs,
         p = DiffEqBase.NullParameters(),
         zero_dummy_der = false,
         kwargs...)
+    inputs isa AbstractVector || (inputs = [inputs])
+    outputs isa AbstractVector || (outputs = [outputs])
     ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs;
         simplify,
         kwargs...)
@@ -1486,7 +1488,7 @@ end
     (; A, B, C, D), simplified_sys = linearize(sys, inputs, outputs;    t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
     (; A, B, C, D)                 = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)
 
-Return a NamedTuple with the matrices of a linear statespace representation
+Linearize `sys` between `inputs` and `outputs`, both vectors of variables. Return a NamedTuple with the matrices of a linear statespace representation
 on the form
 
 ```math
diff --git a/test/linearize.jl b/test/linearize.jl
index 8561d26819..d77793275c 100644
--- a/test/linearize.jl
+++ b/test/linearize.jl
@@ -26,6 +26,13 @@ lsys, ssys = linearize(sys, [r], [r])
 @test lsys.C[] == 0
 @test lsys.D[] == 1
 
+lsys, ssys = linearize(sys, r, r) # Test allow scalars
+
+@test lsys.A[] == -2
+@test lsys.B[] == 1
+@test lsys.C[] == 0
+@test lsys.D[] == 1
+
 ##
 ```