diff --git a/.travis.yml b/.travis.yml
index 367760e..eba13aa 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,8 +4,8 @@ os:
   - linux
   # - osx
 julia:
-  - 0.7
   - 1.0
+  - 1.1
   - nightly
 notifications:
   email: false
diff --git a/Project.toml b/Project.toml
index 33e6e92..9d4c612 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "StaticPolynomials"
 uuid = "62e018b1-6e46-5407-a5a7-97d4fbcae734"
 license = "MIT"
 repo = "https://github.com/JuliaAlgebra/StaticPolynomials.jl.git"
-version = "1.1.1"
+version = "1.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/appveyor.yml b/appveyor.yml
deleted file mode 100644
index fe9c939..0000000
--- a/appveyor.yml
+++ /dev/null
@@ -1,47 +0,0 @@
-environment:
-  matrix:
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe"
-  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
-  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
-
-## uncomment the following lines to allow failures on nightly julia
-## (tests will run but not make your overall status red)
-#matrix:
-#  allow_failures:
-#  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
-#  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
-
-branches:
-  only:
-    - master
-    - /release-.*/
-
-notifications:
-  - provider: Email
-    on_build_success: false
-    on_build_failure: false
-    on_build_status_changed: false
-
-install:
-  - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12"
-# If there's a newer build queued for the same PR, cancel this one
-  - ps: if ($env:APPVEYOR_PULL_REQUEST_NUMBER -and $env:APPVEYOR_BUILD_NUMBER -ne ((Invoke-RestMethod `
-        https://ci.appveyor.com/api/projects/$env:APPVEYOR_ACCOUNT_NAME/$env:APPVEYOR_PROJECT_SLUG/history?recordsNumber=50).builds | `
-        Where-Object pullRequestId -eq $env:APPVEYOR_PULL_REQUEST_NUMBER)[0].buildNumber) { `
-        throw "There are newer queued builds for this pull request, failing early." }
-# Download most recent Julia Windows binary
-  - ps: (new-object net.webclient).DownloadFile(
-        $env:JULIA_URL,
-        "C:\projects\julia-binary.exe")
-# Run installer silently, output to C:\projects\julia
-  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
-
-build_script:
-# Need to convert from shallow to complete for Pkg.clone to work
-  - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
-      Pkg.clone(pwd(), \"StaticPolynomials\"); Pkg.build(\"StaticPolynomials\")"
-
-test_script:
-  - C:\projects\julia\bin\julia -e "Pkg.test(\"StaticPolynomials\")"
diff --git a/docs/Project.toml b/docs/Project.toml
index f4138b3..687aaf2 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,4 +3,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 StaticPolynomials = "62e018b1-6e46-5407-a5a7-97d4fbcae734"
 
 [compat]
-Documenter = "^0.19.6"
+Documenter = "^0.22.0"
diff --git a/docs/make.jl b/docs/make.jl
index cebbf06..92f5d76 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,7 +1,6 @@
 using Documenter, StaticPolynomials
 
 makedocs(
-    format = :html,
     sitename = "StaticPolynomials.jl",
     pages = [
         "Introduction" => "index.md",
@@ -11,9 +10,4 @@ makedocs(
 
 deploydocs(
     repo   = "github.com/JuliaAlgebra/StaticPolynomials.jl.git",
-    target = "build",
-    julia = "0.7",
-    osname = "linux",
-    deps   = nothing,
-    make   = nothing
 )
diff --git a/docs/src/reference.md b/docs/src/reference.md
index 548bffb..e7a9a09 100644
--- a/docs/src/reference.md
+++ b/docs/src/reference.md
@@ -16,6 +16,10 @@ gradient
 gradient!
 evaluate_and_gradient
 evaluate_and_gradient!
+hessian(f::Polynomial, xp...)
+hessian!(U, f::Polynomial, xp...)
+gradient_and_hessian
+gradient_and_hessian!
 differentiate_parameters(::Polynomial, x, p)
 differentiate_parameters!(u, ::Polynomial, x, p)
 ```
@@ -38,6 +42,8 @@ jacobian
 jacobian!
 evaluate_and_jacobian
 evaluate_and_jacobian!
+hessian!(U, F::PolynomialSystem, x::AbstractVector)
+jacobian_and_hessian!
 differentiate_parameters(::PolynomialSystem, x, p)
 differentiate_parameters!(U, ::PolynomialSystem, x, p)
 ```
diff --git a/src/evaluate_codegen.jl b/src/evaluate_codegen.jl
index 44596e4..9092307 100644
--- a/src/evaluate_codegen.jl
+++ b/src/evaluate_codegen.jl
@@ -6,38 +6,38 @@ This assumes that E is in reverse lexicographic order.
 """
 function generate_evaluate(E, ::Type{T}, access_input) where T
     exprs = []
-    last_expr = generate_evaluate!(exprs, E, T, size(E, 1), 1, 1:size(E, 2), access_input)
+    coeffs = map(i -> :(c[$i]), 1:size(E, 2))
+    last_expr = generate_evaluate!(exprs, E, T, size(E, 1), 1, coeffs, access_input)
     Expr(:block, exprs..., last_expr)
 end
 
-function generate_evaluate!(exprs, E, ::Type{T}, nvar, nterm, coeffperm, access_input) where T
+function generate_evaluate!(exprs, E, ::Type{T}, nvar, nterm, coeffs, access_input) where T
     m, n = size(E)
 
     if n == 1
-        return first(monomial_product(T, E[:,1], :(c[$(coeffperm[nterm])]), access_input=access_input))
+        return first(monomial_product(T, E[:,1], coeffs[nterm], access_input=access_input))
     end
 
     if m == 1
-        coeffs = [:(c[$(coeffperm[j])]) for j=nterm:nterm+n-1]
-        return evalpoly(T, E[1,:], coeffs, access_input(nvar))
+        return evalpoly(T, E[1,:], coeffs[nterm:nterm+n-1], access_input(nvar))
     end
 
     # Recursive
-    coeffs = []
+    rec_coeffs = []
     sub_nterm = nterm
     degrees, submatrices = degrees_submatrices(E)
     for E_d in submatrices
-        coeff_subexpr = generate_evaluate!(exprs, E_d, T, nvar - 1, sub_nterm, coeffperm, access_input)
+        coeff_subexpr = generate_evaluate!(exprs, E_d, T, nvar - 1, sub_nterm, coeffs, access_input)
         # If the returned value is just a symbol propagate it
         if isa(coeff_subexpr, Symbol)
-            push!(coeffs, coeff_subexpr)
+            push!(rec_coeffs, coeff_subexpr)
         else
             @gensym c
             push!(exprs, :($c = $coeff_subexpr))
-            push!(coeffs, c)
+            push!(rec_coeffs, c)
         end
         sub_nterm += size(E_d, 2)
     end
 
-    return evalpoly(T, degrees, coeffs, access_input(nvar))
+    return evalpoly(T, degrees, rec_coeffs, access_input(nvar))
 end
diff --git a/src/evaluation.jl b/src/evaluation.jl
index 61146c0..62e5d1e 100644
--- a/src/evaluation.jl
+++ b/src/evaluation.jl
@@ -1,6 +1,7 @@
 export evaluate, gradient, gradient!,
     evaluate_and_gradient, evaluate_and_gradient!,
-    differentiate_parameters, differentiate_parameters!
+    differentiate_parameters, differentiate_parameters!,
+    hessian!, hessian, gradient_and_hessian, gradient_and_hessian!
 
 import Base: @propagate_inbounds
 
@@ -178,6 +179,133 @@ end
     val
 end
 
+
+
+#############
+## Hessian ##
+#############
+
+function _gradient_hessian_impl(f::Type{Polynomial{T, E, P}}; return_grad=true) where {T, E, P}
+    if P == Nothing
+        access_input = x_
+    else
+        n = size(P, 1)
+        access_input(i) = i ≤ n ? :(p[$i]) : :(x[$(i-n)])
+    end
+
+    exps = exponents(E)
+    params = exponents(P)
+
+    n, m = size(E)
+
+    code = []
+    ∇ = [Symbol("∇", i) for i=1:n]
+    vals = [Symbol("val", i) for i=1:n]
+
+
+    # for each variable we simulate the derivative
+    for i in 1:n
+        coeffsᵢ = Union{Symbol, Expr}[]
+        expsᵢ = Vector{eltype(exps)}()
+        paramsᵢ = P === nothing ? nothing : Vector{eltype(exps)}()
+        mᵢ = 0
+
+        # For each term make the symbolic derivative
+        for j in 1:m
+            if exps[i,j] == 0
+                continue
+            end
+            for k in 1:n
+                if i == k
+                    push!(expsᵢ, exps[k, j] - 1)
+                    if exps[i,j] > 1
+                        push!(coeffsᵢ, :($(exps[i,j]) * c[$j]))
+                    else
+                        push!(coeffsᵢ, :(c[$j]))
+                    end
+                else
+                    push!(expsᵢ, exps[k, j])
+                end
+            end
+            if params !== nothing
+                for k in 1:size(params, 1)
+                    push!(paramsᵢ, params[k, j])
+                end
+            end
+
+            mᵢ += 1
+        end
+        if mᵢ == 0
+            zero_tuple = Expr(:tuple, (:(zero($T)) for _=1:n)...)
+            val_grad_codeᵢ = :((zero($T), SVector($zero_tuple)))
+        else
+            Eᵢ = reshape(expsᵢ, (n, mᵢ))
+            Pᵢ = P === nothing ? nothing : reshape(paramsᵢ, (size(paramsᵢ, 1), mᵢ))
+            val_grad_codeᵢ = generate_gradient(Eᵢ, Pᵢ, T, access_input, coeffsᵢ)
+        end
+
+        push!(code, :($(Expr(:tuple, vals[i], ∇[i])) = begin $(val_grad_codeᵢ) end))
+    end
+
+
+    push!(code, :(grad = SVector($(Expr(:tuple, vals...)))))
+    push!(code, :(hess = assemble_matrix(SVector($(Expr(:tuple, (∇[i] for i=1:n)...))))))
+    quote
+        @boundscheck length(x) ≥ $(size(E)[1])
+        c = coefficients(f)
+        $(code...)
+        grad, hess
+    end
+end
+
+@generated function gradient_and_hessian(f::Polynomial{T,E,Nothing}, x::AbstractVector) where {T, E}
+    _gradient_hessian_impl(f)
+end
+@generated function gradient_and_hessian(f::Polynomial, x::AbstractVector, p)
+    _gradient_hessian_impl(f)
+end
+
+@doc """
+     gradient_and_hessian(F::PolynomialSystem, x, [p])
+
+Compute the gradient and Hessian of `f` at `x`.
+""" gradient_and_hessian(F::Polynomial, x)
+
+"""
+    gradient_and_hessian!(u, U, f::Polynomial, x, [p])
+
+Compute the gradient and Hessian of `f` at `x` and store them in `u` and `U`.
+"""
+function gradient_and_hessian!(u, U, f::Polynomial, xp...)
+    grad, hess = gradient_and_hessian(f, xp...)
+    u .= grad
+    U .= hess
+    nothing
+end
+
+"""
+    hessian(f::Polynomial, x, [p])
+
+Compute the hessian of `f` at `x`.
+"""
+function hessian(f::Polynomial, xp...)
+    last(gradient_and_hessian(f, xp...))
+end
+
+"""
+    hessian!(U, f::Polynomial, x, [p])
+
+Compute the hessian of `f` at `x` and store it in `U`.
+"""
+function hessian!(U, f::Polynomial, xp...)
+    U .= hessian(f, xp...)
+    U
+end
+
+
+################
+## Parameters ##
+################
 function _differentiate_parameters_impl(f::Type{Polynomial{T, E, P}}) where {T, E, P}
     @assert P != Nothing
 
diff --git a/src/gradient_codegen.jl b/src/gradient_codegen.jl
index 40c88d0..c17a680 100644
--- a/src/gradient_codegen.jl
+++ b/src/gradient_codegen.jl
@@ -5,32 +5,31 @@ Generate the statements for the evaluation of the polynomial with exponents `E`
 and paramter exponents `P`.
 This assumes that `E` and `P` are in reverse lexicographic order.
 """
-function generate_gradient(E, P, ::Type{T}, access_input) where T
+function generate_gradient(E, P, ::Type{T}, access_input,coeffs=map(i -> :(c[$i]), 1:size(E, 2))) where T
     nvars = size(E,1)
     if P !== nothing
         nvars += size(P, 1)
     end
-    coeffperm = 1:size(E, 2)
 
     exprs = []
     # dvals is a list of tuples (s, iszero) where s is a symbol or expression
     # representing the partial derivative and iszero a boolean flag
     # indicating whether s is zero.
-    val, dvals = partial_derivatives!(exprs, E, P, T, nvars, 1, coeffperm, access_input)
+    val, dvals = partial_derivatives!(exprs, E, P, T, nvars, 1, coeffs, access_input)
     out = :($(val), SVector($(Expr(:tuple, map(first, dvals)...))))
 
     Expr(:block, exprs..., out)
 end
 
-function partial_derivatives!(exprs, E, P, ::Type{T}, nvar, nterm, coeffperm, access_input) where T
+function partial_derivatives!(exprs, E, P, ::Type{T}, nvar, nterm, coeffs, access_input) where T
     m, n = size(E)
 
     # We only have one Term remaining. So we just evaluate the term and compute all
     # partial derivatives
     if n == 1
-        return term(E, P, T, nterm, coeffperm, access_input)
+        return term(E, P, T, nterm, coeffs, access_input)
     elseif m == 1
-        return univariate!(exprs, E, P, T, nvar, nterm, coeffperm, access_input)
+        return univariate!(exprs, E, P, T, nvar, nterm, coeffs, access_input)
     end
 
     degrees, submatrices = degrees_submatrices(E)
@@ -38,7 +37,7 @@ function partial_derivatives!(exprs, E, P, ::Type{T}, nvar, nterm, coeffperm, ac
     dvalues = []
 
     for (k, E_d) in enumerate(submatrices)
-        val, dvals = partial_derivatives!(exprs, E_d, P, T, nvar - 1, nterm, coeffperm, access_input)
+        val, dvals = partial_derivatives!(exprs, E_d, P, T, nvar - 1, nterm, coeffs, access_input)
         push!(values, val)
         push!(dvalues, dvals)
 
@@ -77,9 +76,9 @@ function partial_derivatives!(exprs, E, P, ::Type{T}, nvar, nterm, coeffperm, ac
     val, dvals
 end
 
-function term(E::AbstractMatrix, P, ::Type{T}, nterm, coeffperm, access_input) where T
+function term(E::AbstractMatrix, P, ::Type{T}, nterm, coeffs, access_input) where T
     @assert size(E, 2) == 1
-    c = term_coefficient_with_parameters(T, P, nterm, coeffperm, access_input)
+    c = term_coefficient_with_parameters(T, P, nterm, coeffs, access_input)
     offset = P === nothing ? 0 : size(P, 1)
     monomial_product_with_derivatives(T, E, c, access_input=i -> access_input(i+offset))
 end
@@ -89,24 +88,24 @@ end
 
 Compute the coefficient of the `nterm`-th term with possible parameters.
 """
-function term_coefficient_with_parameters(T, P, nterm, coeffperm, access_input)
-    monomial_product(T, P[:,nterm], :(c[$(coeffperm[nterm])]), access_input=access_input)[1]
+function term_coefficient_with_parameters(T, P, nterm, coeffs, access_input)
+    monomial_product(T, P[:,nterm], coeffs[nterm], access_input=access_input)[1]
 end
-term_coefficient_with_parameters(T, ::Nothing, nterm, coeffperm, access_input) = :(c[$(coeffperm[nterm])])
+term_coefficient_with_parameters(T, ::Nothing, nterm, coeffs, access_input) = coeffs[nterm]
 
 """
     univariate!(exprs, E::AbstractMatrix, P, ::Type{T}, nvar, nterm, access_input)
 
 `E` represents an univariate polynomial.
 """
-function univariate!(exprs, E::AbstractMatrix, P, ::Type{T}, nvar, nterm, coeffperm, access_input) where T
+function univariate!(exprs, E::AbstractMatrix, P, ::Type{T}, nvar, nterm, coeffs, access_input) where T
     @assert size(E, 1) == 1
     n = size(E, 2)
     # If we have parameters it can happen that there are duplicates E
     # For each group of duplicates we have to evaluate the coefficient polynomial
-    coeffs, E_filtered = coefficients_with_parameters!(exprs, E, P, T, nterm, coeffperm, access_input)
+    rec_coeffs, E_filtered = coefficients_with_parameters!(exprs, E, P, T, nterm, coeffs, access_input)
 
-    val, dval = evalpoly_derivative!(exprs, T, E_filtered, coeffs, access_input(nvar))
+    val, dval = evalpoly_derivative!(exprs, T, E_filtered, rec_coeffs, access_input(nvar))
     val, [(dval, n == 1)]
 end
 
@@ -117,13 +116,13 @@ If we have parameters it can happen that there are duplicates in E.
 For each group of duplicates we have to evaluate the coefficient polynomial.
 If there are no parameters this is simply the coefficients vector.
 """
-function coefficients_with_parameters!(exprs, E, ::Nothing, ::Type{T}, nterm, coeffperm, access_input) where T
+function coefficients_with_parameters!(exprs, E, ::Nothing, ::Type{T}, nterm, coeffs, access_input) where T
     n = size(E, 2)
-    [:(c[$(coeffperm[j])]) for j=nterm:nterm+n-1], vec(E)
+    coeffs[nterm:nterm+n-1], vec(E)
 end
-function coefficients_with_parameters!(exprs, E, P, ::Type{T}, nterm, coeffperm, access_input) where T
+function coefficients_with_parameters!(exprs, E, P, ::Type{T}, nterm, coeffs, access_input) where T
     n = size(E, 2)
-    coeffs = []
+    pcoeffs = []
     last = 1
     j = 1
     E_filtered = [E[1,1]]
@@ -141,18 +140,18 @@ function coefficients_with_parameters!(exprs, E, P, ::Type{T}, nterm, coeffperm,
 
         nduplicates = j - last - 1
         if nduplicates == 0 # nothing happened
-            c, _ = monomial_product(T, P[:,nterm_sub], :(c[$(coeffperm[nterm_sub])]), access_input=access_input)
-            push!(coeffs, c)
+            c, _ = monomial_product(T, P[:,nterm_sub], coeffs[nterm_sub], access_input=access_input)
+            push!(pcoeffs, c)
         else
             # we have duplicates, so we consider the rest (the parameters) as an polynomial
             # and pass it to evaluate
             P_sub = @view P[:,nterm_sub:nterm+j-2]
-            c = generate_evaluate!(exprs, P_sub, T, size(P, 1), nterm_sub, coeffperm, access_input)
-            push!(coeffs, c)
+            c = generate_evaluate!(exprs, P_sub, T, size(P, 1), nterm_sub, coeffs, access_input)
+            push!(pcoeffs, c)
         end
         last = j
     end
-    coeffs, E_filtered
+    pcoeffs, E_filtered
 end
 
 function generate_differentiate_parameters(E, P, ::Type{T}, access_input) where T
@@ -164,7 +163,7 @@ function generate_differentiate_parameters(E, P, ::Type{T}, access_input) where
     p = revlexicographic_cols_perm([E; P])
     E = E[:, p]
     P = P[:, p]
-    coeffperm = p
+    coeffs = map(i -> :(c[$i]), p)
     # Since we are only interested in the derivative we can by hand eliminate all
     # terms where no parameter occurs (these are constant terms in our new setting)
     # Since we have a revlexicographic order of the columns, the all 0 column
@@ -190,7 +189,7 @@ function generate_differentiate_parameters(E, P, ::Type{T}, access_input) where
         end
         E = E[:,nconstants+1:end]
         P = P[:,nconstants+1:end]
-        coeffperm = coeffperm[nconstants+1:end]
+        coeffs = coeffs[nconstants+1:end]
     end
     # we can derivate wrt the parameters by changing the roles
     E, P = P, E
@@ -199,7 +198,7 @@ function generate_differentiate_parameters(E, P, ::Type{T}, access_input) where
     # dvals is a list of tuples (s, iszero) where s is a symbol or expression
     # representing the partial derivative and iszero a boolean flag
     # indicating whether s is zero.
-    _, dvals = partial_derivatives!(exprs, E, P, T, nvars, 1, coeffperm, access_input)
+    _, dvals = partial_derivatives!(exprs, E, P, T, nvars, 1, coeffs, access_input)
     quote
         $(exprs...)
         SVector($(Expr(:tuple, map(first, dvals)...)))
diff --git a/src/helpers.jl b/src/helpers.jl
index c975857..17b3168 100644
--- a/src/helpers.jl
+++ b/src/helpers.jl
@@ -71,3 +71,11 @@ function revlexicographic_cols_perm(A::AbstractMatrix; kws...)
     cols = map(i -> (@view A[end:-1:1, i]), inds)
     sortperm(cols; kws...)
 end
+
+@generated function assemble_matrix(vs::SVector{M, SVector{N, T}}) where {T, N, M}
+    quote
+        SMatrix{$M, $N, $T, $(M*N)}(
+            tuple($([:(vs[$i][$j]) for j=1:N for i=1:M]...))
+        )
+    end
+end
diff --git a/src/polynomial_system.jl b/src/polynomial_system.jl
index f3f0c8a..e1c7f4c 100644
--- a/src/polynomial_system.jl
+++ b/src/polynomial_system.jl
@@ -1,6 +1,7 @@
 export PolynomialSystem, npolynomials,
     evaluate!,
     jacobian!, jacobian, evaluate_and_jacobian!, evaluate_and_jacobian,
+    jacobian_and_hessian!,
     differentiate_parameters, differentiate_parameters!
 
 import Base: @propagate_inbounds
@@ -224,15 +225,6 @@ and store its result in `U`.
 Evaluate the Jacobian of the polynomial system `F` at `x` with parameters `p`
 and store its result in `U`.
 """ jacobian!(U, F::PolynomialSystem, x, p)
-
-@generated function assemble_matrix(vs::SVector{M, SVector{N, T}}) where {T, N, M}
-    quote
-        SMatrix{$M, $N, $T, $(M*N)}(
-            tuple($([:(vs[$i][$j]) for j=1:N for i=1:M]...))
-        )
-    end
-end
-
 @generated function _jacobian(F::PolynomialSystem{N}, x...) where {N}
     ∇ = [Symbol("∇", i) for i=1:N]
     quote
@@ -268,7 +260,7 @@ Evaluate the Jacobian of the polynomial system `F` at `x` with parameters `p`.
     quote
         @boundscheck checkbounds(u, 1:$N)
         @boundscheck checkbounds(U, 1:$N, 1:$NVars)
-        
+
         @inbounds begin
             $(map(1:N) do i
                 quote
@@ -334,6 +326,65 @@ Evaluate the system `F` and its Jacobian at `x`.
 Evaluate the system `F` and its Jacobian at `x` with parameters `p`.
 """ evaluate_and_jacobian(F::PolynomialSystem, x, p)
 
+
+############
+# HESSIAN
+###########
+@generated function _hessian!(U::AbstractArray{T,3}, F::PolynomialSystem{N, NVars}, x...) where {T, N, NVars}
+    quote
+        @boundscheck checkbounds(U, 1:$N, 1:$NVars, 1:$NVars)
+        @inbounds begin
+            $(map(1:N) do i
+                quote
+                    H = hessian(F.polys[$i], x...)
+                    @inbounds for k=1:$NVars, j=1:$NVars
+                        U[$i, k, j] = H[k, j]
+                    end
+                end
+            end...)
+        end
+        U
+    end
+end
+
+"""
+     hessian!(U::AbstractArray{T,3}, F::PolynomialSystem, x, [p])
+
+Evaluate the Hessian of the polynomial system `F` at `x` and store its result in `U`.
+"""
+@propagate_inbounds hessian!(U, F::PolynomialSystem, x::AbstractVector) = _hessian!(U, F, x)
+@propagate_inbounds hessian!(U, F::PolynomialSystem, x::AbstractVector, p) = _hessian!(U, F, x, p)
+
+@generated function _jacobian_and_hessian!(u::AbstractMatrix, U::AbstractArray{T,3}, F::PolynomialSystem{N, NVars}, x...) where {T, N, NVars}
+    quote
+        @boundscheck checkbounds(U, 1:$N, 1:$NVars, 1:$NVars)
+        @boundscheck checkbounds(u, 1:$N, 1:$NVars)
+        @inbounds begin
+            $(map(1:N) do i
+                quote
+                    grad, H = gradient_and_hessian(F.polys[$i], x...)
+                    @inbounds for j=1:$NVars
+                        u[$i, j] = grad[j]
+                    end
+                    @inbounds for k=1:$NVars, j=1:$NVars
+                        U[$i, k, j] = H[k, j]
+                    end
+                end
+            end...)
+        end
+        nothing
+    end
+end
+"""
+     jacobian_and_hessian!(u, U::AbstractArray{T,3}, F::PolynomialSystem, x)
+
+Evaluate the Jacobian and Hessian of the polynomial system `F` at `x`
+and store its result in `u` and `U`.
+"""
+@propagate_inbounds jacobian_and_hessian!(u, U, F::PolynomialSystem, x::AbstractVector) = _jacobian_and_hessian!(u, U, F, x)
+@propagate_inbounds jacobian_and_hessian!(u, U, F::PolynomialSystem, x::AbstractVector, p) = _jacobian_and_hessian!(u, U, F, x, p)
+
+
 #####################
 # Parameter Jacobian
 #####################
diff --git a/test/basic_tests.jl b/test/basic_tests.jl
index bc01231..1f62ae7 100644
--- a/test/basic_tests.jl
+++ b/test/basic_tests.jl
@@ -175,6 +175,38 @@
         @test U == [4407 -14297; 1 3]
     end
 
+    @testset "Hessian" begin
+        @polyvar x y z
+
+        w = [2, 3, 5]
+
+        f = x^4*y+x*y*z+y+x+z
+        g = Polynomial(f)
+        grad_g = Polynomial.(MP.differentiate(f, [x, y, z]))
+        hess_w = vcat(map(grad_g) do g_i
+            gradient(g_i, w)'
+        end...)
+
+        @test hessian(g, w) == hess_w
+        @test gradient_and_hessian(g, w) == (gradient(g, w), hess_w)
+
+        U = zeros(Int, 3, 3)
+        @test hessian!(U, g, w) == hess_w
+
+
+        p = [x^4*y+x*y*z+y+x+z, (x^3+y^2)*z]
+        P = PolynomialSystem(p)
+
+        u = zeros(2, 3)
+        U = zeros(2, 3, 3)
+
+        hessian!(U, P, w)
+        @test U[1,:,:] == hess_w
+        jacobian_and_hessian!(u, U, P, w)
+        @test u[1,:] == gradient(g, w)
+        @test U[1,:,:] == hess_w
+    end
+
     @testset "foreach" begin
         @polyvar x y
         g = [Polynomial(x^2+y^2), Polynomial(2x^2+4y^2+3x*y^4+1)]