diff --git a/CMakeLists.txt b/CMakeLists.txt
index b835ec346..c21a78b1a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -261,6 +261,7 @@ set(HEYOKA_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/mascon.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/vsop2013.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/cr3bp.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/model/ffnn.cpp"
     # Math functions.
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepE.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepF.cpp"
diff --git a/include/heyoka/kw.hpp b/include/heyoka/kw.hpp
index c28481287..6cbadce24 100644
--- a/include/heyoka/kw.hpp
+++ b/include/heyoka/kw.hpp
@@ -30,6 +30,13 @@ IGOR_MAKE_NAMED_ARGUMENT(parallel_mode);
 IGOR_MAKE_NAMED_ARGUMENT(prec);
 IGOR_MAKE_NAMED_ARGUMENT(mu);
 
+// kwargs for the ffnn
+IGOR_MAKE_NAMED_ARGUMENT(inputs);
+IGOR_MAKE_NAMED_ARGUMENT(nn_hidden);
+IGOR_MAKE_NAMED_ARGUMENT(n_out);
+IGOR_MAKE_NAMED_ARGUMENT(activations);
+IGOR_MAKE_NAMED_ARGUMENT(nn_wb);
+
 } // namespace kw
 
 HEYOKA_END_NAMESPACE
diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
new file mode 100644
index 000000000..e0141d763
--- /dev/null
+++ b/include/heyoka/model/ffnn.hpp
@@ -0,0 +1,155 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
+//
+// This file is part of the heyoka library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef HEYOKA_MODEL_FFNN_HPP
+#define HEYOKA_MODEL_FFNN_HPP
+
+#include <cstdint>
+#include <functional>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <boost/numeric/conversion/cast.hpp>
+#include <boost/safe_numerics/safe_integer.hpp>
+
+#include <heyoka/config.hpp>
+#include <heyoka/detail/igor.hpp>
+#include <heyoka/detail/type_traits.hpp>
+#include <heyoka/detail/visibility.hpp>
+#include <heyoka/expression.hpp>
+#include <heyoka/kw.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+namespace model
+{
+namespace detail
+{
+
+template <typename... KwArgs>
+auto ffnn_common_opts(const KwArgs &...kw_args)
+{
+    igor::parser p{kw_args...};
+
+    static_assert(!p.has_unnamed_arguments(), "This function accepts only named arguments");
+
+    // Network inputs. Mandatory.
+    // The kw::inputs argument must be a range of values from which
+    // an expression can be constructed.
+    std::vector<expression> inputs;
+    if constexpr (p.has(kw::inputs)) {
+        for (const auto &val : p(kw::inputs)) {
+            inputs.emplace_back(val);
+        }
+    } else {
+        static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                      "The 'inputs' keyword argument is necessary but it was not provided");
+    }
+
+    // Number of hidden neurons per hidden layer. Mandatory.
+    // The kw::nn_hidden argument must be a range containing
+    // integral values.
+    std::vector<std::uint32_t> nn_hidden;
+    if constexpr (p.has(kw::nn_hidden)) {
+        for (const auto &nval : p(kw::nn_hidden)) {
+            nn_hidden.push_back(boost::numeric_cast<std::uint32_t>(nval));
+        }
+    } else {
+        static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                      "The 'nn_hidden' keyword argument is necessary but it was not provided");
+    }
+
+    // Number of network outputs. Mandatory.
+    // The kw::n_out argument must be of integral type.
+    auto n_out = [&p]() {
+        if constexpr (p.has(kw::n_out)) {
+            return boost::numeric_cast<std::uint32_t>(p(kw::n_out));
+        } else {
+            static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                          "The 'n_out' keyword argument is necessary but it was not provided");
+        }
+    }();
+
+    // Network activation functions. Mandatory.
+    // The kw::activations argument must be a range containing values
+    // from which a std::function can be constructed.
+    std::vector<std::function<expression(const expression &)>> activations;
+    if constexpr (p.has(kw::activations)) {
+        for (const auto &f : p(kw::activations)) {
+            activations.emplace_back(f);
+        }
+    } else {
+        static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                      "The 'activations' keyword argument is necessary but it was not provided");
+    }
+
+    // Network weights and biases. Optional, defaults to heyoka parameters.
+    // The kw::nn_wb argument, if present, must be a range of values from which
+    // expressions can be constructed.
+    std::vector<expression> nn_wb;
+    if constexpr (p.has(kw::nn_wb)) {
+        for (const auto &val : p(kw::nn_wb)) {
+            nn_wb.emplace_back(val);
+        }
+    } else {
+        // Safe counterpart to std::uint32_t in order to avoid
+        // overflows when manipulating indices and sizes.
+        using su32 = boost::safe_numerics::safe<std::uint32_t>;
+
+        // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons).
+        auto n_hidden_layers = su32(nn_hidden.size());
+        // Number of neuronal layers (counting input and output).
+        auto n_layers = n_hidden_layers + 2;
+        // Number of inputs.
+        auto n_in = su32(inputs.size());
+        // Number of neurons per neuronal layer.
+        std::vector<su32> n_neurons{n_in};
+        n_neurons.insert(n_neurons.end(), nn_hidden.begin(), nn_hidden.end());
+        n_neurons.insert(n_neurons.end(), n_out);
+
+        // Number of network parameters (wb: weights and biases, w: only weights).
+        su32 n_wb = 0;
+        for (su32 i = 1; i < n_layers; ++i) {
+            n_wb += n_neurons[i - 1] * n_neurons[i];
+            n_wb += n_neurons[i];
+        }
+        nn_wb.resize(n_wb);
+        for (decltype(nn_wb.size()) i = 0; i < nn_wb.size(); ++i) {
+            nn_wb[i] = par[boost::numeric_cast<std::uint32_t>(i)];
+        }
+    }
+
+    return std::tuple{std::move(inputs), std::move(nn_hidden), std::move(n_out), std::move(activations),
+                      std::move(nn_wb)};
+}
+
+// This c++ function returns the symbolic expressions of the `n_out` output neurons in a feed forward neural network,
+// as a function of the `n_in` input expressions.
+//
+// The expression will contain the weights and biases of the neural network flattened into `pars` with the following
+// conventions:
+//
+// from the left to right layer of parameters: [W01, W12,W23, ..., B1,B2,B3,....] where the weight matrices Wij are
+// to be considered as flattened (row first) and so are the bias vectors.
+//
+HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(const std::vector<expression> &, const std::vector<std::uint32_t> &,
+                                                    std::uint32_t,
+                                                    const std::vector<std::function<expression(const expression &)>> &,
+                                                    const std::vector<expression> &);
+} // namespace detail
+
+inline constexpr auto ffnn = [](const auto &...kw_args) -> std::vector<expression> {
+    return std::apply(detail::ffnn_impl, detail::ffnn_common_opts(kw_args...));
+};
+
+} // namespace model
+
+HEYOKA_END_NAMESPACE
+
+#endif
diff --git a/include/heyoka/models.hpp b/include/heyoka/models.hpp
index 12b11a67f..545034687 100644
--- a/include/heyoka/models.hpp
+++ b/include/heyoka/models.hpp
@@ -10,6 +10,7 @@
 #define HEYOKA_MODELS_HPP
 
 #include <heyoka/model/cr3bp.hpp>
+#include <heyoka/model/ffnn.hpp>
 #include <heyoka/model/fixed_centres.hpp>
 #include <heyoka/model/mascon.hpp>
 #include <heyoka/model/nbody.hpp>
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
new file mode 100644
index 000000000..4e5358021
--- /dev/null
+++ b/src/model/ffnn.cpp
@@ -0,0 +1,135 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
+//
+// This file is part of the heyoka library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <algorithm>
+#include <cassert>
+#include <cstdint>
+#include <functional>
+#include <stdexcept>
+#include <vector>
+
+#include <boost/safe_numerics/safe_integer.hpp>
+
+#include <fmt/core.h>
+
+#include <heyoka/config.hpp>
+#include <heyoka/expression.hpp>
+#include <heyoka/model/ffnn.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+namespace model::detail
+{
+
+namespace
+{
+
+// Safe counterpart to std::uint32_t in order to avoid
+// overflows when manipulating indices and sizes.
+using su32 = boost::safe_numerics::safe<std::uint32_t>;
+
+std::vector<expression> compute_layer(su32 layer_id, const std::vector<expression> &inputs,
+                                      const std::vector<su32> &n_neurons,
+                                      const std::function<expression(const expression &)> &activation,
+                                      const std::vector<expression> &nn_wb, su32 n_net_w, su32 &wcounter,
+                                      su32 &bcounter)
+{
+    assert(layer_id > 0u);
+    auto n_neurons_prev_layer = su32(inputs.size());
+    auto n_neurons_curr_layer = n_neurons[layer_id];
+
+    std::vector<expression> retval(static_cast<std::vector<expression>::size_type>(n_neurons_curr_layer), 0_dbl);
+    for (su32 i = 0; i < n_neurons_curr_layer; ++i) {
+        for (su32 j = 0; j < n_neurons_prev_layer; ++j) {
+
+            // Add the weight and update the weight counter.
+            retval[i] += nn_wb[wcounter] * inputs[j];
+            ++wcounter;
+        }
+
+        // Add the bias and update the counter.
+        retval[i] += nn_wb[bcounter + n_net_w];
+        ++bcounter;
+        // Activation function.
+        retval[i] = activation(retval[i]);
+    }
+    return retval;
+}
+
+} // namespace
+
+std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::vector<std::uint32_t> &nn_hidden,
+                                  std::uint32_t n_out,
+                                  const std::vector<std::function<expression(const expression &)>> &activations,
+                                  const std::vector<expression> &nn_wb)
+{
+    // Sanity checks.
+    if (activations.empty()) {
+        throw std::invalid_argument("Cannot create a FFNN with an empty list of activation functions");
+    }
+    if (nn_hidden.size() != activations.size() - 1u) {
+        throw std::invalid_argument(fmt::format(
+            "The number of hidden layers, as detected from the inputs, was {}, while "
+            "the number of activation function supplied was {}. A FFNN needs exactly one more activation function "
+            "than the number of hidden layers.",
+            nn_hidden.size(), activations.size()));
+    }
+    if (in.empty()) {
+        throw std::invalid_argument("The inputs provided to the FFNN is an empty vector.");
+    }
+    if (n_out == 0u) {
+        throw std::invalid_argument("The number of network outputs cannot be zero.");
+    }
+    if (!std::all_of(nn_hidden.begin(), nn_hidden.end(), [](auto item) { return item > 0u; })) {
+        throw std::invalid_argument("The number of neurons for each hidden layer must be greater than zero!");
+    }
+    if (std::any_of(activations.begin(), activations.end(), [](const auto &func) { return !func; })) {
+        throw std::invalid_argument("The list of activation functions cannot contain empty functions");
+    }
+
+    // From now on, always use safe arithmetics to compute/manipulate
+    // indices and sizes.
+    using detail::su32;
+
+    // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons).
+    auto n_hidden_layers = su32(nn_hidden.size());
+    // Number of neuronal layers (counting input and output).
+    auto n_layers = n_hidden_layers + 2;
+    // Number of inputs.
+    auto n_in = su32(in.size());
+    // Number of neurons per neuronal layer.
+    std::vector<su32> n_neurons{n_in};
+    n_neurons.insert(n_neurons.end(), nn_hidden.begin(), nn_hidden.end());
+    n_neurons.insert(n_neurons.end(), n_out);
+    // Number of network parameters (wb: weights and biases, w: only weights).
+    su32 n_net_wb = 0, n_net_w = 0;
+    for (su32 i = 1; i < n_layers; ++i) {
+        n_net_wb += n_neurons[i - 1u] * n_neurons[i];
+        n_net_w += n_neurons[i - 1u] * n_neurons[i];
+        n_net_wb += n_neurons[i];
+    }
+    // Sanity check.
+    if (nn_wb.size() != n_net_wb) {
+        throw std::invalid_argument(fmt::format(
+            "The number of network parameters, detected from its structure to be {}, does not match the size of "
+            "the corresponding expressions: {}.",
+            static_cast<std::uint32_t>(n_net_wb), nn_wb.size()));
+    }
+
+    // Now we build the expressions recursively transvering from layer to layer (L = f(Wx+b))).
+    std::vector<expression> retval = in;
+    su32 wcounter = 0, bcounter = 0;
+    for (su32 i = 1; i < n_layers; ++i) {
+        retval = detail::compute_layer(i, retval, n_neurons, activations[i - 1u], nn_wb, n_net_w, wcounter, bcounter);
+    }
+    return retval;
+}
+
+} // namespace model::detail
+
+HEYOKA_END_NAMESPACE
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 2aef033c8..936c9b9c6 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -142,6 +142,7 @@ ADD_HEYOKA_TESTCASE(model_fixed_centres)
 ADD_HEYOKA_TESTCASE(model_rotating)
 ADD_HEYOKA_TESTCASE(model_mascon)
 ADD_HEYOKA_TESTCASE(model_cr3bp)
+ADD_HEYOKA_TESTCASE(model_ffnn)
 ADD_HEYOKA_TESTCASE(step_callback)
 ADD_HEYOKA_TESTCASE(llvm_state_mem_cache)
 
diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
new file mode 100644
index 000000000..4a5bd8c25
--- /dev/null
+++ b/test/model_ffnn.cpp
@@ -0,0 +1,114 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
+//
+// This file is part of the heyoka library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <cstdint>
+#include <functional>
+#include <limits>
+#include <stdexcept>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <heyoka/config.hpp>
+#include <heyoka/expression.hpp>
+#include <heyoka/kw.hpp>
+#include <heyoka/math.hpp>
+#include <heyoka/model/ffnn.hpp>
+#include <system_error>
+
+#include "catch.hpp"
+
+using namespace heyoka;
+
+TEST_CASE("impl")
+{
+    using Catch::Matchers::Message;
+
+    // A linear layer, just because
+    auto linear = [](expression ret) -> expression { return ret; };
+    // We also define a few symbols
+    auto [x, y, z] = make_vars("x", "y", "z");
+
+    // First, we test malformed cases and their throws.
+    // 1 - number of activations function is wrong
+    REQUIRE_THROWS_AS(model::detail::ffnn_impl({x}, {1}, 2, {heyoka::tanh, heyoka::tanh, linear},
+                                               {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl}),
+                      std::invalid_argument);
+    // 2 - number of inputs is zero
+    REQUIRE_THROWS_AS(
+        model::detail::ffnn_impl({}, {1}, 2, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl}),
+        std::invalid_argument);
+    // 3 - number of outputs is zero
+    REQUIRE_THROWS_AS(model::detail::ffnn_impl({x}, {1}, 0, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 4_dbl}),
+                      std::invalid_argument);
+    // 4 - One of the hidden layers has zero neurons
+    REQUIRE_THROWS_AS(model::detail::ffnn_impl({x}, {1, 0}, 2, {heyoka::tanh, heyoka::tanh, linear}, {1_dbl, 2_dbl}),
+                      std::invalid_argument);
+    // 5 - Wrong number of weights/biases
+    REQUIRE_THROWS_AS(
+        model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 5_dbl, 6_dbl}),
+        std::invalid_argument);
+    // 6 - empty activations list
+    REQUIRE_THROWS_MATCHES(model::detail::ffnn_impl({x}, {1}, 1, {}, {1_dbl, 2_dbl, 3_dbl, 5_dbl, 6_dbl}),
+                           std::invalid_argument,
+                           Message("Cannot create a FFNN with an empty list of activation functions"));
+    // 7 - empty activation functions
+    REQUIRE_THROWS_MATCHES(
+        model::detail::ffnn_impl({x}, {}, 1, {std::function<expression(const expression &)>{}}, {1_dbl, 2_dbl}),
+        std::invalid_argument, Message("The list of activation functions cannot contain empty functions"));
+
+    // We now check some hand coded networks
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {}, 1, {linear}, {1_dbl, 2_dbl});
+        REQUIRE(my_net[0] == expression(2_dbl + x));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {}, 1, {heyoka::tanh}, {1_dbl, 2_dbl});
+        REQUIRE(my_net[0] == expression(heyoka::tanh(2_dbl + x)));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, linear}, {1_dbl, 2_dbl, 3_dbl, 4_dbl});
+        REQUIRE(my_net[0] == expression(4_dbl + (2_dbl * heyoka::tanh(3_dbl + x))));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, heyoka::sin}, {1_dbl, 2_dbl, 3_dbl, 4_dbl});
+        REQUIRE(my_net[0] == expression(heyoka::sin(4_dbl + (2_dbl * heyoka::tanh(3_dbl + x)))));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x, y}, {2}, 1, {heyoka::sin, heyoka::cos},
+                                               {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        REQUIRE(my_net[0] == expression(heyoka::cos(1_dbl + 2_dbl * heyoka::sin(1_dbl + x + y))));
+    }
+
+    // Check overflow detection.
+    REQUIRE_THROWS_AS(
+        model::detail::ffnn_impl({x}, {}, std::numeric_limits<std::uint32_t>::max(), {linear}, {1_dbl, 2_dbl}),
+        std::system_error);
+}
+
+TEST_CASE("igor_iface")
+{
+    auto [x, y, z] = make_vars("x", "y", "z");
+    {
+        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = {2u}, kw::n_out = 1u,
+                                  kw::activations = {heyoka::sin, heyoka::cos},
+                                  kw::nn_wb = {1., 1., 1., 1., 1., 1., 1., 1., 1.});
+        auto vanilla_v = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
+                                                  {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        REQUIRE(igor_v == vanilla_v);
+    }
+    // We test the expected setting for the default weights+biases expressions to par[i].
+    {
+        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = {2u}, kw::n_out = 1u,
+                                  kw::activations = {heyoka::sin, heyoka::cos});
+        auto vanilla_v
+            = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
+                                       {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8]});
+        REQUIRE(igor_v == vanilla_v);
+    }
+}