diff --git a/docs/src/API/electrical.md b/docs/src/API/electrical.md
index cc9dc2b35..788796abf 100644
--- a/docs/src/API/electrical.md
+++ b/docs/src/API/electrical.md
@@ -34,6 +34,8 @@ IdealOpAmp
 Diode
 HeatingDiode
 VariableResistor
+NMOS
+PMOS
 PNP
 NPN
 ```
diff --git a/docs/src/tutorials/MOSFET_calibration.md b/docs/src/tutorials/MOSFET_calibration.md
new file mode 100644
index 000000000..ee7cab059
--- /dev/null
+++ b/docs/src/tutorials/MOSFET_calibration.md
@@ -0,0 +1,74 @@
+# MOSFET I-V Curves 
+In this example, first we'll demonstrate the I-V curves of the NMOS transistor model. 
+First of all, we construct a circuit using the NMOS transistor. We'll need to import ModelingToolkit and the Electrical standard library that holds the transistor models.
+
+```@example NMOS
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t
+using ModelingToolkitStandardLibrary.Electrical
+using ModelingToolkitStandardLibrary.Blocks: Constant
+using OrdinaryDiffEq
+using Plots
+```
+
+Here we just connect the source pin to ground, the drain pin to a voltage source named `Vcc`, and the gate pin to a voltage source named `Vb`.
+```@example NMOS
+@mtkmodel SimpleNMOSCircuit begin
+    @components begin
+        Q1 = NMOS()
+        Vcc = Voltage()
+        Vb = Voltage()
+        ground = Ground()
+
+        Vcc_const = Constant(k=V_cc)
+        Vb_const = Constant(k=V_b)
+    end
+
+    @parameters begin
+        V_cc = 5.0
+        V_b = 3.5
+    end
+    @equations begin
+        #voltage sources
+        connect(Vcc_const.output, Vcc.V)
+        connect(Vb_const.output, Vb.V)
+
+        #ground connections
+        connect(Vcc.n, Vb.n, ground.g, Q1.s)
+
+        #other stuff
+        connect(Vcc.p, Q1.d)
+        connect(Vb.p, Q1.g)
+    end
+end
+
+@mtkbuild sys = SimpleNMOSCircuit(V_cc = 5.0, V_b = 3.5)
+
+prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+sol = solve(prob)
+```
+
+Now to make sure that the transistor model is working like it's supposed to, we can examine the plots of the drain-source voltage vs. the drain current, otherwise knowns as the I-V curve of the transistor. 
+```@example NMOS
+v_cc_list = collect(0.05:0.1:10.0)
+
+I_D_list = []
+I_D_lists = []
+
+for V_b in [1.0, 1.4, 1.8, 2.2, 2.6]
+    I_D_list = []
+    for V_cc in v_cc_list
+        @mtkbuild sys = SimpleNMOSCircuit(V_cc=V_cc, V_b=V_b)
+        prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+        sol = solve(prob)
+        push!(I_D_list, sol[sys.Q1.d.i][1])
+    end
+    push!(I_D_lists, I_D_list)
+end
+
+reduce(hcat, I_D_lists)
+plot(v_cc_list, I_D_lists, title="NMOS IV Curves", label=["V_GS: 1.0 V" "V_GS: 1.4 V" "V_GS: 1.8 V" "V_GS: 2.2 V" "V_GS: 2.6 V"], xlabel = "Drain-Source Voltage (V)", ylabel = "Drain Current (A)")
+```
+
+We can see that we get exactly what we would expect: as the drain-source voltage increases, the drain current increases, until the the transistor gets in to the saturation region of operation. 
+Then the only increase in drain current is due to the channel-length modulation effect. Additionally, we can see that the maximum current reached increases as the gate voltage increases. 
\ No newline at end of file
diff --git a/src/Electrical/Analog/mosfets.jl b/src/Electrical/Analog/mosfets.jl
new file mode 100644
index 000000000..bd2917757
--- /dev/null
+++ b/src/Electrical/Analog/mosfets.jl
@@ -0,0 +1,180 @@
+"""
+    NMOS(;name, V_tn, R_DS, lambda)
+
+Creates an N-type MOSFET transistor 
+
+    # Structural Parameters
+        - `use_transconductance`: If `true` the parameter `k_n` needs to be provided, and is used in the calculation of the current
+        through the transistor. Otherwise, `mu_n`, `C_ox`, `W`, and `L` need to be provided and are used to calculate the transconductance. 
+
+        - `use_channel_length_modulation`: If `true` the channel length modulation effect is taken in to account. In essence this gives  
+        the drain-source current has a small dependency on the drains-source voltage in the saturation region of operation. 
+
+    # Connectors
+        - `d` Drain Pin
+        - `g` Gate Pin
+        - `s` Source Pin
+
+    # Parameters
+        - `mu_n`: Electron mobility
+        - `C_ox`: Oxide capacitance (F/m^2)
+        - `W`: Channel width (m)
+        - `L`: Channel length
+        - `k_n`: MOSFET transconductance parameter 
+
+Based on the MOSFET models in (Sedra, A. S., Smith, K. C., Carusone, T. C., & Gaudet, V. C. (2021). Microelectronic circuits (8th ed.). Oxford University Press.)
+"""
+@mtkmodel NMOS begin
+    @variables begin
+        V_GS(t)
+        V_DS(t)
+        V_OV(t)
+    end
+
+    @components begin
+        d = Pin()
+        g = Pin()
+        s = Pin()
+    end
+
+    @parameters begin
+        V_tn = 0.8, [description = "Threshold voltage (V)"]
+        R_DS = 1e7, [description = "Drain to source resistance (Ω)"]
+
+        if use_channel_length_modulation
+            lambda = 1 / 25, [description = "Channel length modulation coefficient (V^(-1))"]
+        end
+
+        if !use_transconductance
+            mu_n, [description = "Electron mobility"]
+            C_ox, [description = "Oxide capacitance (F/m^2)"]
+            W, [description = "Channel width (m)"]
+            L, [description = "Channel length (m)"]
+        else
+            k_n = 20e-3, [description = "MOSFET transconductance parameter"]
+        end
+    end
+
+    @structural_parameters begin
+        use_transconductance = true
+        use_channel_length_modulation = true
+    end
+
+    begin
+        if !use_transconductance
+            k_n = mu_n * C_ox * (W / L)
+        end
+    end
+
+    @equations begin
+        V_DS ~ ifelse(d.v < s.v, s.v - d.v, d.v - s.v)
+        V_GS ~ g.v - ifelse(d.v < s.v, d.v, s.v)
+        V_OV ~ V_GS - V_tn
+
+        d.i ~ ifelse(d.v < s.v, -1, 1) * ifelse(V_GS < V_tn,
+            0 + V_DS / R_DS,
+            ifelse(V_DS < V_OV,
+                ifelse(use_channel_length_modulation,
+                    k_n * (1 + lambda * V_DS) * (V_OV - V_DS / 2) * V_DS + V_DS / R_DS,
+                    k_n * (V_OV - V_DS / 2) * V_DS + V_DS / R_DS),
+                ifelse(use_channel_length_modulation,
+                    ((k_n * V_OV^2) / 2) * (1 + lambda * V_DS) + V_DS / R_DS,
+                    (k_n * V_OV^2) / 2 + V_DS / R_DS)
+            )
+        )
+
+        g.i ~ 0
+        s.i ~ -d.i
+    end
+end
+
+
+
+"""
+    PMOS(;name, V_tp, R_DS, lambda)
+
+Creates an N-type MOSFET transistor 
+
+    # Structural Parameters
+        - `use_transconductance`: If `true` the parameter `k_p` needs to be provided, and is used in the calculation of the current
+        through the transistor. Otherwise, `mu_n`, `C_ox`, `W`, and `L` need to be provided and are used to calculate the transconductance. 
+
+        - `use_channel_length_modulation`: If `true` the channel length modulation effect is taken in to account. In essence this gives  
+        the drain-source current has a small dependency on the drains-source voltage in the saturation region of operation. 
+
+    # Connectors
+        - `d` Drain Pin
+        - `g` Gate Pin
+        - `s` Source Pin
+
+    # Parameters
+        - `mu_p`: Electron mobility
+        - `C_ox`: Oxide capacitance (F/m^2)
+        - `W`: Channel width (m)
+        - `L`: Channel length
+        - `k_p`: MOSFET transconductance parameter 
+
+Based on the MOSFET models in (Sedra, A. S., Smith, K. C., Carusone, T. C., & Gaudet, V. C. (2021). Microelectronic circuits (8th ed.). Oxford University Press.)
+"""
+@mtkmodel PMOS begin
+    @variables begin
+        V_GS(t)
+        V_DS(t)
+    end
+
+    @components begin
+        d = Pin()
+        g = Pin()
+        s = Pin()
+    end
+
+    @parameters begin
+        V_tp = -1.5, [description = "Threshold voltage (V)"]
+        R_DS = 1e7, [description = "Drain-source resistance (Ω)"]
+
+        if use_channel_length_modulation
+            lambda = 1 / 25, [description = "Channel length modulation coefficient (V^(-1))"]
+        end
+
+        if !use_transconductance
+            mu_p, [description = "Hole mobility"]
+            C_ox, [description = "Oxide capacitance (F/m^2)"]
+            W, [description = "Channel width (m)"]
+            L, [description = "Channel length (m)"]
+        else
+            k_p = 20e-3
+        end
+    end
+
+    @structural_parameters begin
+        use_transconductance = true
+        use_channel_length_modulation = true
+    end
+
+    begin
+        if !use_transconductance
+            k_p = mu_p * C_ox * (W / L)
+        end
+    end
+
+    @equations begin
+        V_DS ~ ifelse(d.v > s.v, s.v - d.v, d.v - s.v)
+        V_GS ~ g.v - ifelse(d.v > s.v, d.v, s.v)
+
+        d.i ~ -ifelse(d.v > s.v, -1.0, 1.0) * ifelse(V_GS > V_tp,
+            0.0 + V_DS / R_DS,
+            ifelse(V_DS > (V_GS - V_tp),
+                ifelse(use_channel_length_modulation,
+                    k_p * (1 + lambda * V_DS) * ((V_GS - V_tp) - V_DS / 2) * V_DS +
+                    V_DS / R_DS,
+                    k_p * ((V_GS - V_tp) - V_DS / 2) * V_DS + V_DS / R_DS),
+                ifelse(use_channel_length_modulation,
+                    ((k_p * (V_GS - V_tp)^2) / 2) * (1 + lambda * V_DS) + V_DS / R_DS,
+                    (k_p * (V_GS - V_tp)^2) / 2 + V_DS / R_DS)
+            )
+        )
+
+        g.i ~ 0
+        s.i ~ -d.i
+    end
+end
\ No newline at end of file
diff --git a/src/Electrical/Electrical.jl b/src/Electrical/Electrical.jl
index baac5138e..712e4b172 100644
--- a/src/Electrical/Electrical.jl
+++ b/src/Electrical/Electrical.jl
@@ -24,6 +24,9 @@ include("Analog/sensors.jl")
 export Voltage, Current
 include("Analog/sources.jl")
 
+export NMOS, PMOS
+include("Analog/mosfets.jl")
+
 export NPN, PNP
 include("Analog/transistors.jl")
 
diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl
index b3f96515b..38b258445 100644
--- a/test/Electrical/analog.jl
+++ b/test/Electrical/analog.jl
@@ -541,6 +541,213 @@ end
     # savefig(plt, "rc_circuit_test_variable_resistor")
 end
 
+@testset "NMOS Transistor" begin
+
+    @mtkmodel SimpleNMOSCircuit begin
+        @components begin
+            Q1 = NMOS()
+            Vcc = Voltage()
+            Vb = Voltage()
+            ground = Ground()
+
+            Vcc_const = Constant(k=V_cc)
+            Vb_const = Constant(k=V_b)
+        end
+
+        @parameters begin
+            V_cc = 5.0
+            V_b = 3.5
+        end
+        @equations begin
+            #voltage sources
+            connect(Vcc_const.output, Vcc.V)
+            connect(Vb_const.output, Vb.V)
+
+            #ground connections
+            connect(Vcc.n, Vb.n, ground.g, Q1.s)
+
+            #other stuff
+            connect(Vcc.p, Q1.d)
+            connect(Vb.p, Q1.g)
+        end
+    end
+
+    @mtkbuild sys = SimpleNMOSCircuit(V_cc = 5.0, V_b = 3.5)
+
+    prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+    sol = solve(prob)
+    @test sol[sys.Q1.d.i][1] > 0.0
+    @test sol[sys.Q1.s.i][1] < 0.0
+    @test sol[sys.Q1.g.i][1] == 0.0
+    @test sol[sys.Q1.d.v][1] == 5.0
+    @test sol[sys.Q1.s.v] < sol[sys.Q1.d.v]
+
+    # test device symmetry
+    @mtkmodel FlippedNMOSCircuit begin
+        @components begin
+            Q1 = NMOS()
+            Vcc = Voltage()
+            Vb = Voltage()
+            ground = Ground()
+
+            Vcc_const = Constant(k=V_cc)
+            Vb_const = Constant(k=V_b)
+        end
+
+        @parameters begin
+            V_cc = 5.0
+            V_b = 3.5
+        end
+        @equations begin
+            #voltage sources
+            connect(Vcc_const.output, Vcc.V)
+            connect(Vb_const.output, Vb.V)
+
+            #ground connections
+            connect(Vcc.n, Vb.n, ground.g, Q1.d)
+
+            #other stuff
+            connect(Vcc.p, Q1.s)
+            connect(Vb.p, Q1.g)
+        end
+    end
+
+    @mtkbuild flipped_sys = FlippedNMOSCircuit(V_cc=5.0, V_b=3.5)
+
+    flipped_prob = ODEProblem(flipped_sys, Pair[], (0.0, 10.0))
+    flipped_sol = solve(flipped_prob)
+    @test flipped_sol[flipped_sys.Q1.d.i][1] < 0
+    @test flipped_sol[flipped_sys.Q1.s.i][1] > 0 
+    @test flipped_sol[flipped_sys.Q1.s.v] > flipped_sol[flipped_sys.Q1.d.v]
+
+    # channel length modulation
+    @mtkmodel SimpleNMOSCircuitChannel begin
+            @components begin
+                Q1 = NMOS(use_channel_length_modulation = false)
+                Vcc = Voltage()
+                Vb = Voltage()
+                ground = Ground()
+
+                Vcc_const = Constant(k=V_cc)
+                Vb_const = Constant(k=V_b)
+            end
+
+            @parameters begin
+                V_cc = 5.0
+                V_b = 3.5
+            end
+            @equations begin
+                #voltage sources
+                connect(Vcc_const.output, Vcc.V)
+                connect(Vb_const.output, Vb.V)
+
+                #ground connections
+                connect(Vcc.n, Vb.n, ground.g, Q1.s)
+
+                #other stuff
+                connect(Vcc.p, Q1.d)
+                connect(Vb.p, Q1.g)
+            end
+    end
+
+        @mtkbuild sys = SimpleNMOSCircuitChannel(V_cc = 5.0, V_b = 3.5)
+
+        prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+        sol = solve(prob)
+        @test sol[sys.Q1.d.i][1] > 0.0
+        @test sol[sys.Q1.s.i][1] < 0.0
+end
+
+
+@testset "PMOS Transistor" begin
+
+    @mtkmodel SimplePMOSCircuit begin
+        @components begin
+            Q1 = PMOS(use_channel_length_modulation = false)
+            Vs = Voltage()
+            Vb = Voltage()
+            Vd = Voltage()
+            ground = Ground()
+
+            Vs_const = Constant(k=V_s)
+            Vb_const = Constant(k=V_b)
+            Vd_const = Constant(k=V_d)
+        end
+
+        @parameters begin
+            V_s = 5.0
+            V_b = 3.5
+            V_d = 0.0
+        end
+        @equations begin
+            #voltage sources
+            connect(Vs_const.output, Vs.V)
+            connect(Vb_const.output, Vb.V)
+            connect(Vd_const.output, Vd.V )
+
+            #ground connections
+            connect(Vs.n, Vb.n, ground.g, Vd.n)
+
+            connect(Vd.p, Q1.d)
+            #other stuff
+            connect(Vs.p, Q1.s)
+            connect(Vb.p, Q1.g)
+        end
+    end
+
+    @mtkbuild sys = SimplePMOSCircuit(V_s=5.0, V_b=2.5, V_d = 3)
+
+    prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+    sol = solve(prob)
+
+    @test sol[sys.Q1.d.i][1] < 0.0
+    @test sol[sys.Q1.s.i][1] > 0.0
+
+    # device symmetry
+    @mtkmodel FlippedPMOSCircuit begin
+        @components begin
+            Q1 = PMOS()
+            Vs = Voltage()
+            Vb = Voltage()
+            Vd = Voltage()
+            ground = Ground()
+
+            Vs_const = Constant(k=V_s)
+            Vb_const = Constant(k=V_b)
+            Vd_const = Constant(k=V_d)
+        end
+
+        @parameters begin
+            V_s = 5.0
+            V_b = 3.5
+            V_d = 0.0
+        end
+        @equations begin
+            #voltage sources
+            connect(Vs_const.output, Vs.V)
+            connect(Vb_const.output, Vb.V)
+            connect(Vd_const.output, Vd.V )
+
+            #ground connections
+            connect(Vs.n, Vb.n, ground.g, Vd.n)
+
+            connect(Vd.p, Q1.s)
+            #other stuff
+            connect(Vs.p, Q1.d)
+            connect(Vb.p, Q1.g)
+        end
+    end
+
+    @mtkbuild flipped_sys = FlippedPMOSCircuit(V_s=5.0, V_b=2.5, V_d = 3)
+
+    flipped_prob = ODEProblem(flipped_sys, Pair[], (0.0, 10.0))
+    flipped_sol = solve(flipped_prob)
+
+    flipped_sol[flipped_sys.Q1.d.i][1] > 0.0
+    flipped_sol[flipped_sys.Q1.s.i][1] < 0.0
+
+
+end
 
 @testset "NPN Tests" begin
 
@@ -710,4 +917,3 @@ end
         sol[sys.Q1.s.i][15],
         0.0,
         atol = 1e-16)
-end
\ No newline at end of file