diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 2d2a04a..9e0d8f8 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -50,6 +50,9 @@ jobs:
     - name: Run test
       run: fpm test
 
+    - name: Run examples
+      run: fpm run --example --all
+
   Docs:
     runs-on: ${{ matrix.os }}
     strategy:
diff --git a/examples/example_hybrd.f90 b/examples/example_hybrd.f90
new file mode 100644
index 0000000..b6bdd81
--- /dev/null
+++ b/examples/example_hybrd.f90
@@ -0,0 +1,103 @@
+!> The problem is to determine the values of x(1), x(2), ..., x(9)
+!>  which solve the system of tridiagonal equations
+!>     (3-2*x(1))*x(1)           -2*x(2)                   = -1
+!>             -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
+!>                                 -x(8) + (3-2*x(9))*x(9) = -1
+program example_hybrd
+
+    use minpack_module, only: hybrd, enorm, dpmpar
+    implicit none
+    integer j, n, maxfev, ml, mu, mode, nprint, info, nfev, ldfjac, lr, nwrite
+    double precision xtol, epsfcn, factor, fnorm
+    double precision x(9), fvec(9), diag(9), fjac(9, 9), r(45), qtf(9), &
+        wa1(9), wa2(9), wa3(9), wa4(9)
+
+    !> Logical output unit is assumed to be number 6.
+    data nwrite/6/
+
+    n = 9
+
+    !> The following starting values provide a rough solution.
+    do j = 1, 9
+        x(j) = -1.0d0
+    end do
+
+    ldfjac = 9
+    lr = 45
+
+    !> Set xtol to the square root of the machine precision.
+    !>  unless high precision solutions are required,
+    !>  this is the recommended setting.
+    xtol = dsqrt(dpmpar(1))
+
+    maxfev = 2000
+    ml = 1
+    mu = 1
+    epsfcn = 0.0d0
+    mode = 2
+    do j = 1, 9
+        diag(j) = 1.0d0
+    end do
+    factor = 1.0d2
+    nprint = 0
+
+    call hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, &
+               mode, factor, nprint, info, nfev, fjac, ldfjac, &
+               r, lr, qtf, wa1, wa2, wa3, wa4)
+    fnorm = enorm(n, fvec)
+    write (nwrite, 1000) fnorm, nfev, info, (x(j), j=1, n)
+
+1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
+           5x, "NUMBER OF FUNCTION EVALUATIONS", i10// &
+           5x, "EXIT PARAMETER", 16x, i10// &
+           5x, "FINAL APPROXIMATE SOLUTION"//(5x, 3d15.7))
+
+    !> Results obtained with different compilers or machines
+    !>  may be slightly different.
+    !>
+    !>> FINAL L2 NORM OF THE RESIDUALS  0.1192636D-07
+    !>>
+    !>> NUMBER OF FUNCTION EVALUATIONS        14
+    !>>
+    !>> EXIT PARAMETER                         1
+    !>>
+    !>> FINAL APPROXIMATE SOLUTION
+    !>>
+    !>>  -0.5706545D+00 -0.6816283D+00 -0.7017325D+00
+    !>>  -0.7042129D+00 -0.7013690D+00 -0.6918656D+00
+    !>>  -0.6657920D+00 -0.5960342D+00 -0.4164121D+00
+
+contains
+
+    !> Subroutine fcn for hybrd example.
+    subroutine fcn(n, x, fvec, iflag)
+
+        implicit none
+        integer, intent(in) :: n
+        integer, intent(inout) :: iflag
+        double precision, intent(in) :: x(n)
+        double precision, intent(out) :: fvec(n)
+
+        integer k
+        double precision one, temp, temp1, temp2, three, two, zero
+        data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
+
+        if (iflag /= 0) go to 5
+
+        !! Insert print statements here when nprint is positive.
+
+        return
+5       continue
+        do k = 1, n
+            temp = (three - two*x(k))*x(k)
+            temp1 = zero
+            if (k /= 1) temp1 = x(k - 1)
+            temp2 = zero
+            if (k /= n) temp2 = x(k + 1)
+            fvec(k) = temp - temp1 - two*temp2 + one
+        end do
+        return
+
+    end subroutine fcn
+
+end program example_hybrd
diff --git a/examples/example_hybrd1.f90 b/examples/example_hybrd1.f90
new file mode 100644
index 0000000..21ae354
--- /dev/null
+++ b/examples/example_hybrd1.f90
@@ -0,0 +1,80 @@
+!> the problem is to determine the values of x(1), x(2), ..., x(9),
+!>  which solve the system of tridiagonal equations.
+!>
+!>  (3-2*x(1))*x(1)           -2*x(2)                   = -1
+!>          -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
+!>                              -x(8) + (3-2*x(9))*x(9) = -1
+program example_hybrd1
+
+    use minpack_module, only: hybrd1, dpmpar, enorm
+    implicit none
+    integer j, n, info, lwa, nwrite
+    double precision tol, fnorm
+    double precision x(9), fvec(9), wa(180)
+
+    !> Logical output unit is assumed to be number 6.
+    data nwrite/6/
+
+    n = 9
+
+    !> The following starting values provide a rough solution.
+    do j = 1, 9
+        x(j) = -1.d0
+    end do
+
+    lwa = 180
+
+    !> Set tol to the square root of the machine precision.
+    !>  unless high precision solutions are required,
+    !>  this is the recommended setting.
+    tol = dsqrt(dpmpar(1))
+
+    call hybrd1(fcn, n, x, fvec, tol, info, wa, lwa)
+    fnorm = enorm(n, fvec)
+    write (nwrite, 1000) fnorm, info, (x(j), j=1, n)
+
+1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
+           5x, "EXIT PARAMETER", 16x, i10// &
+           5x, "FINAL APPROXIMATE SOLUTION"// &
+           (5x, 3d15.7))
+
+    !> Results obtained with different compilers or machines
+    !>  may be slightly different.
+    !>
+    !>> FINAL L2 NORM OF THE RESIDUALS  0.1192636D-07
+    !>>
+    !>> EXIT PARAMETER                         1
+    !>>
+    !>> FINAL APPROXIMATE SOLUTION
+    !>>
+    !>>  -0.5706545D+00 -0.6816283D+00 -0.7017325D+00
+    !>>  -0.7042129D+00 -0.7013690D+00 -0.6918656D+00
+    !>>  -0.6657920D+00 -0.5960342D+00 -0.4164121D+00
+
+contains
+
+    !> Subroutine fcn for hybrd1 example.
+    subroutine fcn(n, x, fvec, iflag)
+
+        implicit none
+        integer, intent(in) :: n
+        integer, intent(inout) :: iflag
+        double precision, intent(in) :: x(n)
+        double precision, intent(out) :: fvec(n)
+
+        integer k
+        double precision one, temp, temp1, temp2, three, two, zero
+        data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
+
+        do k = 1, n
+            temp = (three - two*x(k))*x(k)
+            temp1 = zero
+            if (k /= 1) temp1 = x(k - 1)
+            temp2 = zero
+            if (k /= n) temp2 = x(k + 1)
+            fvec(k) = temp - temp1 - two*temp2 + one
+        end do
+
+    end subroutine fcn
+
+end program example_hybrd1
diff --git a/examples/example_lmder1.f90 b/examples/example_lmder1.f90
new file mode 100644
index 0000000..8b235bb
--- /dev/null
+++ b/examples/example_lmder1.f90
@@ -0,0 +1,93 @@
+module testmod_der1
+implicit none
+private
+public fcn, dp
+
+integer, parameter :: dp=kind(0d0)
+
+contains
+
+subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag)
+integer, intent(in) :: m, n, ldfjac
+integer, intent(inout) :: iflag
+real(dp), intent(in) :: x(n)
+real(dp), intent(inout) :: fvec(m), fjac(ldfjac, n)
+
+integer :: i
+real(dp) :: tmp1, tmp2, tmp3, tmp4, y(15)
+y = [1.4D-1, 1.8D-1, 2.2D-1, 2.5D-1, 2.9D-1, 3.2D-1, 3.5D-1, 3.9D-1, &
+        3.7D-1, 5.8D-1, 7.3D-1, 9.6D-1, 1.34D0, 2.1D0, 4.39D0]
+
+if (iflag == 1) then
+    do i = 1, 15
+        tmp1 = i
+        tmp2 = 16 - i
+        tmp3 = tmp1
+        if (i .gt. 8) tmp3 = tmp2
+        fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
+    end do
+else
+    do i = 1, 15
+        tmp1 = i
+        tmp2 = 16 - i
+        tmp3 = tmp1
+        if (i .gt. 8) tmp3 = tmp2
+        tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
+        fjac(i,1) = -1.D0
+        fjac(i,2) = tmp1*tmp2/tmp4
+        fjac(i,3) = tmp1*tmp3/tmp4
+    end do
+end if
+end subroutine
+
+end module
+
+
+program example_lmder1
+use minpack_module, only: enorm, lmder1, chkder
+use testmod_der1, only: dp, fcn
+implicit none
+
+integer :: info
+real(dp) :: tol, x(3), fvec(15), fjac(size(fvec), size(x))
+integer :: ipvt(size(x))
+real(dp), allocatable :: wa(:)
+
+! The following starting values provide a rough fit.
+x = [1._dp, 1._dp, 1._dp]
+
+call check_deriv()
+
+! Set tol to the square root of the machine precision. Unless high precision
+! solutions are required, this is the recommended setting.
+tol = sqrt(epsilon(1._dp))
+
+allocate(wa(5*size(x) + size(fvec)))
+call lmder1(fcn, size(fvec), size(x), x, fvec, fjac, size(fjac, 1), tol, &
+    info, ipvt, wa, size(wa))
+print 1000, enorm(size(fvec), fvec), info, x
+1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
+            5x, 'EXIT PARAMETER', 16x, i10              // &
+            5x, 'FINAL APPROXIMATE SOLUTION'            // &
+            5x, 3d15.7)
+
+contains
+
+subroutine check_deriv()
+integer :: iflag
+real(dp) :: xp(size(x)), fvecp(size(fvec)), err(size(fvec))
+call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, &
+    1, err)
+iflag = 1
+call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag)
+iflag = 2
+call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag)
+iflag = 1
+call fcn(size(fvec), size(x), xp, fvecp, fjac, size(fjac, 1), iflag)
+call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, &
+    2, err)
+print *, "Derivatives check (1.0 is correct, 0.0 is incorrect):"
+print *, err
+end subroutine
+
+end program
diff --git a/examples/example_lmdif1.f90 b/examples/example_lmdif1.f90
new file mode 100644
index 0000000..8b508b0
--- /dev/null
+++ b/examples/example_lmdif1.f90
@@ -0,0 +1,61 @@
+module testmod_dif1
+implicit none
+private
+public fcn, dp
+
+integer, parameter :: dp=kind(0d0)
+
+contains
+
+subroutine fcn(m, n, x, fvec, iflag)
+integer, intent(in) :: m, n
+integer, intent(inout) :: iflag
+real(dp), intent(in) :: x(n)
+real(dp), intent(out) :: fvec(m)
+
+integer :: i
+real(dp) :: tmp1, tmp2, tmp3, y(15)
+! Suppress compiler warning:
+y(1) = iflag
+y = [1.4D-1, 1.8D-1, 2.2D-1, 2.5D-1, 2.9D-1, 3.2D-1, 3.5D-1, 3.9D-1, &
+        3.7D-1, 5.8D-1, 7.3D-1, 9.6D-1, 1.34D0, 2.1D0, 4.39D0]
+
+do i = 1, 15
+    tmp1 = i
+    tmp2 = 16 - i
+    tmp3 = tmp1
+    if (i .gt. 8) tmp3 = tmp2
+    fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
+end do
+end subroutine
+
+end module
+
+
+program example_lmdif1
+use minpack_module, only: enorm, lmdif1
+use testmod_dif1, only: dp, fcn
+implicit none
+
+integer :: info, m, n
+real(dp) :: tol, x(3), fvec(15)
+integer :: iwa(size(x))
+real(dp), allocatable :: wa(:)
+
+! The following starting values provide a rough fit.
+x = [1._dp, 1._dp, 1._dp]
+
+! Set tol to the square root of the machine precision. Unless high precision
+! solutions are required, this is the recommended setting.
+tol = sqrt(epsilon(1._dp))
+
+m = size(fvec)
+n = size(x)
+allocate(wa(m*n + 5*n + m))
+call lmdif1(fcn, size(fvec), size(x), x, fvec, tol, info, iwa, wa, size(wa))
+print 1000, enorm(size(fvec), fvec), info, x
+1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
+            5x, 'EXIT PARAMETER', 16x, i10              // &
+            5x, 'FINAL APPROXIMATE SOLUTION'            // &
+            5x, 3d15.7)
+end program
diff --git a/examples/example_primes.f90 b/examples/example_primes.f90
new file mode 100644
index 0000000..c014bd3
--- /dev/null
+++ b/examples/example_primes.f90
@@ -0,0 +1,96 @@
+module types
+implicit none
+private
+public dp
+
+integer, parameter :: dp=kind(0d0)
+
+end module
+
+module find_fit_module
+
+! This module contains a general function find_fit() for a nonlinear least
+! squares fitting. The function can fit any nonlinear expression to any data.
+
+use minpack_module, only: lmdif1
+use types, only: dp
+implicit none
+private
+public find_fit
+
+contains
+
+subroutine find_fit(data_x, data_y, expr, pars)
+! Fits the (data_x, data_y) arrays with the function expr(x, pars).
+! The user can provide any nonlinear function 'expr' depending on any number of
+! parameters 'pars' and it must return the evaluated expression on the
+! array 'x'. The arrays 'data_x' and 'data_y' must have the same
+! length.
+real(dp), intent(in) :: data_x(:), data_y(:)
+interface
+    function expr(x, pars) result(y)
+    use types, only: dp
+    implicit none
+    real(dp), intent(in) :: x(:), pars(:)
+    real(dp) :: y(size(x))
+    end function
+end interface
+real(dp), intent(inout) :: pars(:)
+
+real(dp) :: tol, fvec(size(data_x))
+integer :: iwa(size(pars)), info, m, n
+real(dp), allocatable :: wa(:)
+
+tol = sqrt(epsilon(1._dp))
+m = size(fvec)
+n = size(pars)
+allocate(wa(m*n + 5*n + m))
+call lmdif1(fcn, m, n, pars, fvec, tol, info, iwa, wa, size(wa))
+if (info /= 1) stop "failed to converge"
+
+contains
+
+subroutine fcn(m, n, x, fvec, iflag)
+integer, intent(in) :: m, n
+integer, intent(inout) :: iflag
+real(dp), intent(in) :: x(n)
+real(dp), intent(out) :: fvec(m)
+! Suppress compiler warning:
+fvec(1) = iflag
+fvec = data_y - expr(data_x, x)
+end subroutine
+
+end subroutine
+
+end module
+
+
+program example_primes
+
+! Find a nonlinear fit of the form a*x*log(b + c*x) to a list of primes.
+
+use find_fit_module, only: find_fit
+use types, only: dp
+implicit none
+
+real(dp) :: pars(3)
+real(dp), parameter :: y(*) = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, &
+    37, 41, 43, 47, 53, 59, 61, 67, 71]
+integer :: i
+pars = [1._dp, 1._dp, 1._dp]
+call find_fit([(real(i, dp), i=1,size(y))], y, expression, pars)
+print *, pars
+
+contains
+
+function expression(x, pars) result(y)
+real(dp), intent(in) :: x(:), pars(:)
+real(dp) :: y(size(x))
+real(dp) :: a, b, c
+a = pars(1)
+b = pars(2)
+c = pars(3)
+y = a*x*log(b + c*x)
+end function
+
+end program
diff --git a/fpm.toml b/fpm.toml
index e79beb8..6d478fc 100644
--- a/fpm.toml
+++ b/fpm.toml
@@ -46,3 +46,28 @@ main = "test_lmdif.f90"
 name = "test_lmstr"
 source-dir = "test"
 main = "test_lmstr.f90"
+
+[[example]]
+name = "example_hybrd"
+source-dir = "examples"
+main = "example_hybrd.f90"
+
+[[example]]
+name = "example_hybrd1"
+source-dir = "examples"
+main = "example_hybrd1.f90"
+
+[[example]]
+name = "example_lmder1"
+source-dir = "examples"
+main = "example_lmder1.f90"
+
+[[example]]
+name = "example_lmdif1"
+source-dir = "examples"
+main = "example_lmdif1.f90"
+
+[[example]]
+name = "example_primes"
+source-dir = "examples"
+main = "example_primes.f90"