diff --git a/Project.toml b/Project.toml
index 2ee0b813..52ef20aa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,7 @@ uuid = "ca8b7239-ccd3-5cce-807f-2072f3f0d108"
 version = "0.9.3-DEV"
 
 [deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
@@ -13,6 +14,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
+Compat = "3.7"
 IterTools = "1"
 OrderedCollections = "1"
 ProgressMeter = "0.6, 0.7, 0.8, 0.9, 1"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 3629e36f..00da65c6 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -1,15 +1,25 @@
 [[ACME]]
-deps = ["IterTools", "LinearAlgebra", "Markdown", "OrderedCollections", "ProgressMeter", "SparseArrays", "StaticArrays", "Statistics"]
+deps = ["Compat", "IterTools", "LinearAlgebra", "Markdown", "OrderedCollections", "ProgressMeter", "SparseArrays", "StaticArrays", "Statistics"]
 path = ".."
 uuid = "ca8b7239-ccd3-5cce-807f-2072f3f0d108"
 
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "2e23d71ad695ec28ca58ddd44869f07afa33cc76"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.7.0"
+
 [[Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
 [[Distributed]]
 deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -101,6 +111,10 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
diff --git a/src/ACME.jl b/src/ACME.jl
index 112e25b9..cb701d89 100644
--- a/src/ACME.jl
+++ b/src/ACME.jl
@@ -1,10 +1,11 @@
-# Copyright 2015, 2016, 2017, 2018, 2019 Martin Holters
+# Copyright 2015, 2016, 2017, 2018, 2019, 2020 Martin Holters
 # See accompanying license file.
 
 module ACME
 
 export DiscreteModel, run!, steadystate, steadystate!, linearize, ModelRunner
 
+using Compat: evalpoly
 using SparseArrays: SparseMatrixCSC, blockdiag, dropzeros!, findnz,
     nonzeros, sparse, spzeros
 using LinearAlgebra: BLAS, I, axpy!, lu, rmul!
diff --git a/src/elements.jl b/src/elements.jl
index df9b42a4..5ecbd306 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -1,4 +1,4 @@
-# Copyright 2015, 2016, 2017, 2018, 2019 Martin Holters
+# Copyright 2015, 2016, 2017, 2018, 2019, 2020 Martin Holters
 # See accompanying license file.
 
 export resistor, potentiometer, capacitor, inductor, transformer,
@@ -394,23 +394,31 @@ Pins: `base`, `emitter`, `collector`
 end
 
 @doc doc"""
-    mosfet(typ; vt=0.7, α=2e-5)
+    mosfet(typ; vt=0.7, α=2e-5, λ=0)
 
 Creates a MOSFET transistor with the simple model
 
 $i_D=\begin{cases}
   0 & \text{if } v_{GS} \le v_T \\
   \alpha \cdot (v_{GS} - v_T - \tfrac{1}{2}v_{DS})\cdot v_{DS}
+  \cdot (1 + \lambda v_{DS})
   & \text{if } v_{DS} \le v_{GS} - v_T \cap v_{GS} > v_T \\
-  \frac{\alpha}{2} \cdot (v_{GS} - v_T)^2 & \text{otherwise.}
+  \frac{\alpha}{2} \cdot (v_{GS} - v_T)^2 \cdot (1 + \lambda v_{DS})
+  & \text{otherwise.}
 \end{cases}$
 
 The `typ` parameter chooses between NMOS (`:n`) and PMOS (`:p`). The threshold
 voltage `vt` is given in Volt, `α` (in A/V²) in a constant depending on the
-physics and dimensions of the device.
+physics and dimensions of the device, and `λ` (in Ω⁻¹) controls the channel
+length modulation.
+
+Optionally, it is possible to specify tuples of coefficients for `vt` and `α`.
+These will be used as polynomials in $v_{GS}$ to determine $v_T$ and $\alpha$,
+respectively. E.g. with `vt=(0.7, 0.1, 0.02)`, the $v_{GS}$-dpendent threshold
+voltage $v_T = 0.7 + 0.1\cdot v_{GS} + 0.02\cdot v_{GS}^2$ will be used.
 
 Pins: `gate`, `source`, `drain`
-""" function mosfet(typ; vt=0.7, α=2e-5)
+""" function mosfet(typ; vt=0.7, α=2e-5, λ=0)
     if typ == :n
         polarity = 1
     elseif typ == :p
@@ -418,25 +426,43 @@ Pins: `gate`, `source`, `drain`
     else
         throw(ArgumentError("Unknown mosfet type $(typ), must be :n or :p"))
     end
-    return Element(mv=[-1 0; 0 -1; 0 0; 0 0],
-        mi=[0 0; 0 0; 0 -1; 1 0],
-        mq=polarity*[1 0 0; 0 1 0; 0 0 1; 0 0 0],
-        u0=polarity*[-vt; 0; 0; 0],
-        ports=[:gate => :source, :drain => :source],
-        nonlinear_eq = @inline function (q)
-            vg, vds, id=q # vg = vgs-vt
-            if vg <= 0
-                res = @SVector [-id]
-                J = @SMatrix [0.0 0.0 -1.0]
-            elseif vds <= vg
-                res = @SVector [α * (vg-0.5*vds)*vds - id]
-                J = @SMatrix [α*vds α*(vg-vds) -1.0]
-            else # 0 < vg < vds
-                res = @SVector [(α/2) * vg^2 - id]
-                J = @SMatrix [α*vg 0.0 -1.0]
-            end
-            return (res, J)
-        end)
+    vt = (vt...,)
+    α = (α...,)
+    dvt = vt[2:end] .* (1:length(vt)-1...,)
+    dα = α[2:end] .* (1:length(α)-1...,)
+    let polarity = polarity, α = α, vt = vt
+        return Element(mv=[-1 0; 0 -1; 0 0; 0 0],
+            mi=[0 0; 0 0; 0 -1; 1 0],
+            mq=polarity*[1 0 0; 0 1 0; 0 0 1; 0 0 0],
+            ports=[:gate => :source, :drain => :source],
+            nonlinear_eq = @inline function (q)
+                vgs, vds, id = q
+                α´ = evalpoly(polarity*vgs, α)
+                if !isempty(dα)
+                    dα´_dvgs = evalpoly(polarity*vgs, dα)
+                else
+                    dα´_dvgs = 0
+                end
+                vt´ = evalpoly(polarity*vgs, vt)
+                if !isempty(dvt)
+                    dvt´_dvgs = evalpoly(polarity*vgs, dvt)
+                else
+                    dvt´_dvgs = 0
+                end
+                λ´ = vds ≥ 0 ? λ : zero(λ)
+                if vgs <= vt´
+                    res = @SVector [-id]
+                    J = @SMatrix [0.0 0.0 -1.0]
+                elseif vds <= vgs - vt´ # && vgs > vt´
+                    res = @SVector [α´ * (vgs-vt´-0.5*vds)*vds*(1+λ´*vds) - id]
+                    J = @SMatrix [α´*(1-dvt´_dvgs)*vds*(1+λ´*vds) + dα´_dvgs * (vgs-vt´-0.5*vds)*vds*(1+λ´*vds)  α´*(vgs-vt´+vds*(2*λ´*(vgs-vt´-0.75*vds)-1))  -1.0]
+                else # 0 < vgs - vt´ < vds
+                    res = @SVector [(α´/2) * (vgs-vt´)^2*(1+λ´*vds) - id]
+                    J = @SMatrix [α´*(vgs-vt´)*(1-dvt´_dvgs)*(1+λ´*vds) + dα´_dvgs/2 * (vgs-vt´)^2*(1+λ´*vds) λ´*α´/2*(vgs-vt´)^2 -1.0]
+                end
+                return (res, J)
+            end)
+    end
 end
 
 @doc doc"""
diff --git a/test/runtests.jl b/test/runtests.jl
index 1c269ad4..b5748896 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,9 +1,10 @@
-# Copyright 2015, 2016, 2017, 2018, 2019 Martin Holters
+# Copyright 2015, 2016, 2017, 2018, 2019, 2020 Martin Holters
 # See accompanying license file.
 
 include("checklic.jl")
 
 using ACME
+using Compat: evalpoly, only
 using Test: @test, @test_broken, @test_logs, @test_throws, @testset
 using FFTW: rfft
 using ProgressMeter
@@ -518,6 +519,28 @@ end
         y = run!(model, pol*[0 1 2 2 2; 5 5 0.5 1 1.5])
         @test y == pol*[0 0 1e-4*(1-0.5/2)*0.5 1e-4*(1-1/2)*1 1e-4/2*1^2]
     end
+    for (typ, pol) in ((:n, 1), (:p, -1)), α in (1e-4, (0.0205, -0.0017)),
+        vt in (1, (1.2078, 0.3238), (-1.2454, -0.199, -0.0483))
+        circ = @circuit begin
+            vgs = voltagesource(), [-] == gnd
+            vds = voltagesource(), [-] == gnd
+            J = mosfet(typ, vt=vt, α=α, λ=0.05), [gate] == vgs[+], [drain] == vds[+]
+            out = currentprobe(), [+] == J[source], [-] == gnd
+        end
+        model = DiscreteModel(circ, 1);
+        for vgs in range(0, stop=5, length=10), vds in range(0, stop=5, length=10)
+            y = only(run!(model, pol*hcat([vgs; vds])))
+            α´ = evalpoly(pol*vgs, (α...,))
+            vt´ = evalpoly(pol*vgs, (vt...,))
+            if vgs ≤ vt´
+                @test y == 0
+            elseif vds ≤ vgs - vt´
+                @test y ≈ pol * α´ * (vgs - vt´ - vds / 2) * vds * (1 + 0.05 * vds)
+            else
+                @test y ≈ pol * α´ / 2 * (vgs - vt´)^2 * (1 + 0.05 * vds)
+            end
+        end
+    end
 end
 
 @testset "op amp" begin