diff --git a/applications/test/ODE/Make/options b/applications/test/ODE/Make/options
index 6a97b686ac6283055ce29c5069b1af9dc855c050..1476c5b0e030400637d98b32cb7a745bc99c69a4 100644
--- a/applications/test/ODE/Make/options
+++ b/applications/test/ODE/Make/options
@@ -1,2 +1,3 @@
-EXE_INC = -I$(LIB_SRC)/ODE/lnInclude
+EXE_INC = /* -DFoam_no_besselFunc */ -I$(LIB_SRC)/ODE/lnInclude
+
 EXE_LIBS = -lODE
diff --git a/src/OpenFOAM/primitives/Scalar/Scalar.H b/src/OpenFOAM/primitives/Scalar/Scalar.H
index c72cb3563c3b4cdf3ae67b488caa89586c192513..833f73fe933a7942282a09446965c606ee5566c4 100644
--- a/src/OpenFOAM/primitives/Scalar/Scalar.H
+++ b/src/OpenFOAM/primitives/Scalar/Scalar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -184,6 +184,8 @@ besselFunc(j0)
 besselFunc(j1)
 besselFunc(y0)
 besselFunc(y1)
+besselFunc2(jn)
+besselFunc2(yn)
 
 
 inline Scalar& setComponent(Scalar& s, const direction)
diff --git a/src/OpenFOAM/primitives/Scalar/doubleScalar/doubleScalar.H b/src/OpenFOAM/primitives/Scalar/doubleScalar/doubleScalar.H
index baa0cb15b3bf83ed27f64e54aae6352dc8054ca3..2f7ed678bb2440ddc724dbfe3a7e06c45b42fdd9 100644
--- a/src/OpenFOAM/primitives/Scalar/doubleScalar/doubleScalar.H
+++ b/src/OpenFOAM/primitives/Scalar/doubleScalar/doubleScalar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2015 OpenFOAM Foundation
@@ -85,31 +85,42 @@ inline Scalar atan2(const Scalar y, const Scalar x)
     return ::atan2(y, x);
 }
 
-inline Scalar jn(const int n, const Scalar s)
-{
-    return ::jn(n, s);
-}
-
-inline Scalar yn(const int n, const Scalar s)
-{
-    return ::yn(n, s);
-}
 
 // Normal (double-precision) transcendental functions
-#define transFunc(func)            \
-inline Scalar func(const Scalar s) \
-{                                  \
-    return ::func(s);              \
-}
-
-
-// Double-precision bessel functions
-#define besselFunc(func)           \
-inline Scalar func(const Scalar s) \
-{                                  \
-    return ::func(s);              \
+#define transFunc(func)                                         \
+inline Scalar func(const Scalar s)                              \
+{                                                               \
+    return ::func(s);                                           \
 }
 
+// Normal (double-precision) bessel functions.
+// May not be available on all systems
+#ifdef Foam_no_besselFunc
+    // Not available
+    #define besselFunc(func)                                    \
+    inline Scalar func(const Scalar s)                          \
+    {                                                           \
+        std::cerr<< "No '" << #func << "' function\n";          \
+        return 0;                                               \
+    }
+    #define besselFunc2(func)                                   \
+    inline Scalar func(const int n, const Scalar s)             \
+    {                                                           \
+        std::cerr<< "No '" << #func << "' function\n";          \
+        return 0;                                               \
+    }
+#else
+    #define besselFunc(func)                                    \
+    inline Scalar func(const Scalar s)                          \
+    {                                                           \
+        return ::func(s);                                       \
+    }
+    #define besselFunc2(func)                                   \
+    inline Scalar func(const int n, const Scalar s)             \
+    {                                                           \
+        return ::func(n, s);                                    \
+    }
+#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -127,6 +138,7 @@ inline Scalar func(const Scalar s) \
 #undef ScalarRead
 #undef transFunc
 #undef besselFunc
+#undef besselFunc2
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/primitives/Scalar/floatScalar/floatScalar.H b/src/OpenFOAM/primitives/Scalar/floatScalar/floatScalar.H
index 9d0cbf7d577518497843166ae7f3db37e392f834..c1030a4ef5e9222c0500b491ef6d8bc42de3f451 100644
--- a/src/OpenFOAM/primitives/Scalar/floatScalar/floatScalar.H
+++ b/src/OpenFOAM/primitives/Scalar/floatScalar/floatScalar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2015 OpenFOAM Foundation
@@ -86,50 +86,54 @@ inline Scalar atan2(const Scalar y, const Scalar x)
 }
 
 // Single-precision transcendental functions (with 'f' appended to the name)
-#define transFunc(func)            \
-inline Scalar func(const Scalar s) \
-{                                  \
-    return ::func##f(s);           \
+#define transFunc(func)                                         \
+inline Scalar func(const Scalar s)                              \
+{                                                               \
+    return ::func##f(s);                                        \
 }
 
 
-#ifdef darwin
-// Single-precision bessel functions. (No float version for darwin).
-#define besselFunc(func)           \
-inline Scalar func(const Scalar s) \
-{                                  \
-    return ::func(s);              \
-}
-
-inline Scalar jn(const int n, const Scalar s)
-{
-    return Scalar(::jn(n, double(s)));
-}
-
-inline Scalar yn(const int n, const Scalar s)
-{
-    return Scalar(::yn(n, double(s)));
-}
+// Single-precision bessel functions.
+// May not be available on all systems
+#ifdef Foam_no_besselFunc
+    // Not available
+    #define besselFunc(func)                                    \
+    inline Scalar func(const Scalar s)                          \
+    {                                                           \
+        std::cerr<< "No '" << #func << "' function\n";          \
+        return 0;                                               \
+    }
+    #define besselFunc2(func)                                   \
+    inline Scalar func(const int n, const Scalar s)             \
+    {                                                           \
+        std::cerr<< "No '" << #func << "' function\n";          \
+        return 0;                                               \
+    }
+#elif defined(darwin)
+    // No float version for darwin - use a cast.
+    #define besselFunc(func)                                    \
+    inline Scalar func(const Scalar s)                          \
+    {                                                           \
+        return ::func(s);                                       \
+    }
+    #define besselFunc2(func)                                   \
+    inline Scalar func(const int n, const Scalar s)             \
+    {                                                           \
+        return Scalar(::func(n, double(s)));                    \
+    }
 
 #else
-
-// Single-precision bessel functions (with 'f' appended to the name)
-#define besselFunc(func)           \
-inline Scalar func(const Scalar s) \
-{                                  \
-    return ::func##f(s);           \
-}
-
-inline Scalar jn(const int n, const Scalar s)
-{
-    return ::jnf(n, s);
-}
-
-inline Scalar yn(const int n, const Scalar s)
-{
-    return ::ynf(n, s);
-}
-
+    // With 'f' (float) appended to the name
+    #define besselFunc(func)                                    \
+    inline Scalar func(const Scalar s)                          \
+    {                                                           \
+        return ::func##f(s);                                    \
+    }
+    #define besselFunc2(func)                                   \
+    inline Scalar func(const int n, const Scalar s)             \
+    {                                                           \
+        return ::func##f(n, s);                                 \
+    }
 #endif
 
 
@@ -149,6 +153,7 @@ inline Scalar yn(const int n, const Scalar s)
 #undef ScalarRead
 #undef transFunc
 #undef besselFunc
+#undef besselFunc2
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //