diff --git a/Project.toml b/Project.toml
index 40e2092712..4919de1aca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.3.0"
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
diff --git a/docs/.gitignore b/docs/.gitignore
new file mode 100644
index 0000000000..a303fff203
--- /dev/null
+++ b/docs/.gitignore
@@ -0,0 +1,2 @@
+build/
+site/
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000000..5961dbf7c4
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,11 @@
+using Documenter, ModelingToolkit, DiffEqBase
+
+
+makedocs(
+    sitename="ModelingToolkit.jl",
+    modules=[ModelingToolkit],
+    pages=[
+        "Home" => "index.md",
+        "api.md",
+    ]
+)
diff --git a/docs/src/api.md b/docs/src/api.md
new file mode 100644
index 0000000000..9646a4ed39
--- /dev/null
+++ b/docs/src/api.md
@@ -0,0 +1,55 @@
+# API
+
+
+## Intermediate representation
+
+### Types
+
+```@docs
+Expression
+Variable
+ModelingToolkit.Constant
+Operation
+Equation
+Differential
+```
+
+### Functions
+
+```@docs
+Base.get(c::ModelingToolkit.Constant)
+Base.:~(::Expression, ::Expression)
+expand_derivatives
+ModelingToolkit.derivative
+```
+
+### Macros
+
+```@docs
+@parameters
+@variables
+@derivatives
+@register
+```
+
+## Systems
+
+### Types
+
+```@docs
+ModelingToolkit.AbstractSystem
+ODESystem
+NonlinearSystem
+```
+
+### Functions
+
+```@docs
+independent_variables
+dependent_variables
+parameters
+calculate_jacobian
+generate_jacobian
+generate_function
+DiffEqBase.ODEFunction(sys::ODESystem, dvs, ps; version::FunctionVersion = ArrayFunction)
+```
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000000..52570560a5
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,12 @@
+# ModelingToolkit.jl Documentation
+
+
+ModelingToolkit.jl is an intermediate representation (IR) of computational graphs
+for scientific computing problems. Its purpose is to be a common target for
+modeling DSLs in order to allow for a common platform for model inspection and
+transformation. It uses a tagged variable IR in order to allow specification of
+complex models and allow for transformations of models. It has ways to plug into
+its function registration and derivative system so that way it can interact
+nicely with user-defined routines. Together, this is an abstract form of a
+scientific model that is easy for humans to generate but also easy for programs
+to manipulate.
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index e5cf37f202..7a6b6e16a2 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -12,20 +12,70 @@ using StaticArrays, LinearAlgebra
 using MacroTools
 import MacroTools: splitdef, combinedef
 
+using DocStringExtensions
+
+"""
+$(TYPEDEF)
+
+Base type for a symbolic expression.
+"""
 abstract type Expression <: Number end
+
+"""
+$(TYPEDEF)
+
+TODO
+"""
 abstract type AbstractSystem end
 
 Base.promote_rule(::Type{<:Number},::Type{<:Expression}) = Expression
 Base.zero(::Type{<:Expression}) = Constant(0)
 Base.one(::Type{<:Expression}) = Constant(1)
 
+"""
+$(TYPEDSIGNATURES)
+
+Calculate the jacobian matrix of a system.
+
+Returns a matrix of [`Expression`](@ref) instances. The result from the first
+call will be cached in the system object.
+"""
 function calculate_jacobian end
+
+"""
+$(TYPEDSIGNATURES)
+
+Generate a function to calculate the Jacobian of the system.
+"""
 function generate_jacobian end
+
+"""
+$(TYPEDSIGNATURES)
+
+Generate a function to evaluate the system's equations.
+"""
 function generate_function end
 
+"""
+$(TYPEDSIGNATURES)
+
+Get the set of independent variables for the given system.
+"""
 function independent_variables end
-function dependent_variables   end
-function parameters            end
+
+"""
+$(TYPEDSIGNATURES)
+
+Get the set of dependent variables for the given system.
+"""
+function dependent_variables end
+
+"""
+$(TYPEDSIGNATURES)
+
+Get the set of parameters variables for the given system.
+"""
+function parameters end
 
 @enum FunctionVersion ArrayFunction=1 SArrayFunction=2
 
diff --git a/src/differentials.jl b/src/differentials.jl
index e0ccf5d9b5..bac00808c0 100644
--- a/src/differentials.jl
+++ b/src/differentials.jl
@@ -1,7 +1,30 @@
 export Differential, expand_derivatives, @derivatives
 
 
+"""
+$(TYPEDEF)
+
+Represents a differential operator.
+
+# Fields
+$(FIELDS)
+
+# Examples
+
+```jldoctest
+julia> using ModelingToolkit
+
+julia> @variables x y;
+
+julia> D = Differential(x)
+(D'~x())
+
+julia> D(y)  # Differentiate y wrt. x
+(D'~x())(y())
+```
+"""
 struct Differential <: Function
+    """The variable or expression to differentiate with respect to."""
     x::Expression
 end
 (D::Differential)(x) = Operation(D, Expression[x])
@@ -11,6 +34,11 @@ Base.convert(::Type{Expr}, D::Differential) = D
 
 Base.:(==)(D1::Differential, D2::Differential) = isequal(D1.x, D2.x)
 
+"""
+$(SIGNATURES)
+
+TODO
+"""
 function expand_derivatives(O::Operation)
     @. O.args = expand_derivatives(O.args)
 
@@ -32,6 +60,40 @@ end
 expand_derivatives(x) = x
 
 # Don't specialize on the function here
+"""
+$(SIGNATURES)
+
+Calculate the derivative of the op `O` with respect to its argument with index
+`idx`.
+
+# Examples
+
+```jldoctest label1
+julia> using ModelingToolkit
+
+julia> @variables x y;
+
+julia> ModelingToolkit.derivative(sin(x), 1)
+cos(x())
+```
+
+Note that the function does not recurse into the operation's arguments, i.e. the
+chain rule is not applied:
+
+```jldoctest label1
+julia> myop = sin(x) * y^2
+sin(x()) * y() ^ 2
+
+julia> typeof(myop.op)  # Op is multiplication function
+typeof(*)
+
+julia> ModelingToolkit.derivative(myop, 1)  # wrt. sin(x)
+y() ^ 2
+
+julia> ModelingToolkit.derivative(myop, 2)  # wrt. y^2
+sin(x())
+```
+"""
 derivative(O::Operation, idx) = derivative(O.op, (O.args...,), Val(idx))
 
 # Pre-defined derivatives
@@ -76,10 +138,33 @@ function _differential_macro(x)
     ex
 end
 
+"""
+$(SIGNATURES)
+
+Define one or more differentials.
+
+# Examples
+
+```jldoctest
+julia> using ModelingToolkit
+
+julia> @variables x y z;
+
+julia> @derivatives Dx'~x Dy'~y  # Create differentials wrt. x and y
+((D'~x()), (D'~y()))
+
+julia> Dx(z)  # Differentiate z wrt. x
+(D'~x())(z())
+
+julia> Dy(z)  # Differentiate z wrt. y
+(D'~y())(z())
+```
+"""
 macro derivatives(x...)
     esc(_differential_macro(x))
 end
 
+
 function calculate_jacobian(eqs, dvs)
     Expression[Differential(dv)(eq) for eq ∈ eqs, dv ∈ dvs]
 end
diff --git a/src/equations.jl b/src/equations.jl
index ffaa222021..db50c51acf 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -1,12 +1,42 @@
 export Equation
 
 
+"""
+$(TYPEDEF)
+
+An equality relationship between two expressions.
+
+# Fields
+$(FIELDS)
+"""
 struct Equation
+    """The expression on the left hand side of the equation."""
     lhs::Expression
+    """The expression on the right hand side of the equation."""
     rhs::Expression
 end
 Base.:(==)(a::Equation, b::Equation) = isequal((a.lhs, a.rhs), (b.lhs, b.rhs))
 
+"""
+$(TYPEDSIGNATURES)
+
+Create an [`Equation`](@ref) out of two [`Expression`](@ref) instances, or an
+`Expression` and a `Number`.
+
+# Examples
+
+```jldoctest
+julia> using ModelingToolkit
+
+julia> @variables x y;
+
+julia> x ~ y
+Equation(x(), y())
+
+julia> x - y ~ 0
+Equation(x() - y(), ModelingToolkit.Constant(0))
+```
+"""
 Base.:~(lhs::Expression, rhs::Expression) = Equation(lhs, rhs)
 Base.:~(lhs::Expression, rhs::Number    ) = Equation(lhs, rhs)
 Base.:~(lhs::Number    , rhs::Expression) = Equation(lhs, rhs)
diff --git a/src/function_registration.jl b/src/function_registration.jl
index 61a2e8a15a..de71d4bb15 100644
--- a/src/function_registration.jl
+++ b/src/function_registration.jl
@@ -1,4 +1,9 @@
 # Register functions and handle literals
+"""
+$(SIGNATURES)
+
+TODO
+"""
 macro register(sig)
     splitsig = splitdef(:($sig = nothing))
     name = splitsig[:name]
diff --git a/src/operations.jl b/src/operations.jl
index eb36c1c979..45afdc78b5 100644
--- a/src/operations.jl
+++ b/src/operations.jl
@@ -1,5 +1,47 @@
+"""
+$(TYPEDEF)
+
+An expression representing the application of a function to symbolic arguments.
+
+# Fields
+$(FIELDS)
+
+# Examples
+
+Operations can be built by application of most built-in mathematical functions
+to other [`Expression`](@ref) instances:
+
+```jldoctest
+julia> using ModelingToolkit
+
+julia> @variables x y;
+
+julia> op1 = sin(x)
+sin(x())
+
+julia> typeof(op1.op)
+typeof(sin)
+
+julia> op1.args
+1-element Array{Expression,1}:
+ x()
+
+julia> op2 = x + y
+x() + y()
+
+julia> typeof(op2.op)
+typeof(+)
+
+julia> op2.args
+2-element Array{Expression,1}:
+ x()
+ y()
+```
+"""
 struct Operation <: Expression
+    """The function to be applied."""
     op::Function
+    """The arguments the function is applied to."""
     args::Vector{Expression}
 end
 
diff --git a/src/systems/diffeqs/diffeqsystem.jl b/src/systems/diffeqs/diffeqsystem.jl
index abb2bc79ab..470b6b7c8e 100644
--- a/src/systems/diffeqs/diffeqsystem.jl
+++ b/src/systems/diffeqs/diffeqsystem.jl
@@ -31,11 +31,43 @@ function to_diffeq(eq::Equation)
 end
 Base.:(==)(a::DiffEq, b::DiffEq) = isequal((a.x, a.n, a.rhs), (b.x, b.n, b.rhs))
 
+"""
+$(TYPEDEF)
+
+A system of ordinary differential equations.
+
+# Fields
+* `eqs` - The ODEs defining the system.
+
+# Examples
+
+```
+using ModelingToolkit
+
+@parameters t σ ρ β
+@variables x(t) y(t) z(t)
+@derivatives D'~t
+
+eqs = [D(x) ~ σ*(y-x),
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
+
+de = ODESystem(eqs)
+```
+"""
 struct ODESystem <: AbstractSystem
+    """The ODEs defining the system."""
     eqs::Vector{DiffEq}
+    """Independent variable."""
     iv::Variable
+    """Dependent (state) variables."""
     dvs::Vector{Variable}
+    """Parameter variables."""
     ps::Vector{Variable}
+    """
+    Jacobian matrix. Note: this field will not be defined until
+    [`calculate_jacobian`](@ref) is called on the system.
+    """
     jac::RefValue{Matrix{Expression}}
 end
 
@@ -156,6 +188,13 @@ function generate_factorized_W(sys::ODESystem, simplify=true; version::FunctionV
     return (Wfact_func, Wfact_t_func)
 end
 
+"""
+$(SIGNATURES)
+
+Create an `ODEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and `ps`
+are used to set the order of the dependent variable and parameter vectors,
+respectively.
+"""
 function DiffEqBase.ODEFunction(sys::ODESystem, dvs, ps; version::FunctionVersion = ArrayFunction)
     expr = generate_function(sys, dvs, ps; version = version)
     if version === ArrayFunction
diff --git a/src/systems/nonlinear/nonlinear_system.jl b/src/systems/nonlinear/nonlinear_system.jl
index fda2626fe9..10580226e4 100644
--- a/src/systems/nonlinear/nonlinear_system.jl
+++ b/src/systems/nonlinear/nonlinear_system.jl
@@ -10,9 +10,32 @@ function Base.convert(::Type{NLEq}, eq::Equation)
 end
 Base.:(==)(a::NLEq, b::NLEq) = a.rhs == b.rhs
 
+"""
+$(TYPEDEF)
+
+A nonlinear system of equations.
+
+# Fields
+* `eqs` - Vector of equations defining the system.
+
+# Examples
+
+```
+@variables x y z
+@parameters σ ρ β
+
+eqs = [0 ~ σ*(y-x),
+       0 ~ x*(ρ-z)-y,
+       0 ~ x*y - β*z]
+ns = NonlinearSystem(eqs, [x,y,z])
+```
+"""
 struct NonlinearSystem <: AbstractSystem
+    """Vector of equations defining the system."""
     eqs::Vector{NLEq}
+    """Unknown variables."""
     vs::Vector{Expression}
+    """Parameters."""
     ps::Vector{Variable}
     function NonlinearSystem(eqs, vs)
         rhss = [eq.rhs for eq ∈ eqs]
diff --git a/src/variables.jl b/src/variables.jl
index 96c8cc6a07..5d729db47c 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -1,8 +1,21 @@
 export Variable, @variables, @parameters
 
 
+"""
+$(TYPEDEF)
+
+A named variable which represents a numerical value. The variable's value may
+be known (parameters, independent variables) or unknown (dependent variables).
+
+# Fields
+$(FIELDS)
+"""
 struct Variable <: Function
+    """The variable's unique name."""
     name::Symbol
+    """
+    Whether the variable's value is known.
+    """
     known::Bool
     Variable(name; known = false) = new(name, known)
 end
@@ -17,9 +30,23 @@ function Base.show(io::IO, ::MIME"text/plain", x::Variable)
 end
 
 
+"""
+$(TYPEDEF)
+
+An expression which wraps a constant numerical value.
+
+The value of the constant can be extracted with [`Base.get`](@ref).
+"""
 struct Constant <: Expression
+    """The constant's numerical value"""
     value::Number
 end
+
+"""
+$(TYPEDSIGNATURES)
+
+Get the value of a [`ModelingToolkit.Constant`](@ref).
+"""
 Base.get(c::Constant) = c.value
 
 Base.iszero(ex::Expression) = isa(ex, Constant) && iszero(ex.value)
@@ -69,9 +96,22 @@ function _parse_vars(macroname, known, x)
     push!(ex.args, build_expr(:tuple, var_names))
     return ex
 end
+
+
+"""
+$(SIGNATURES)
+
+Define one or more unknown variables.
+"""
 macro variables(xs...)
     esc(_parse_vars(:variables, false, xs))
 end
+
+"""
+$(SIGNATURES)
+
+Define one or more known variables.
+"""
 macro parameters(xs...)
     esc(_parse_vars(:parameters, true, xs))
 end