diff --git a/embed/README.md b/embed/README.md
new file mode 100644
index 0000000..c824fe1
--- /dev/null
+++ b/embed/README.md
@@ -0,0 +1,21 @@
+# Embedding Python in C
+
+Goal: to make the existing RADAE Python implementation a C library
+
+Approaches:
+
+1. Cython: docs suggest the top level must be a Python program
+2. [Embedding](https://docs.python.org/3/extending/embedding.html): C main() call Python functions. Currently testing this approach 
+3. libtorch: worth exploring
+
+# Packages
+
+`sudo apt install python3.10-dev`
+
+# Build and Run docs.python.org demo
+
+```
+gcc embed1.c -o embed1 $(python3.10-config --clags) $(python3.10-config --ldflags --embed)
+PYTHONPATH="." ./embed1 mult multiply 2 2
+```
+
diff --git a/embed/embed1.c b/embed/embed1.c
new file mode 100644
index 0000000..ecc7fb6
--- /dev/null
+++ b/embed/embed1.c
@@ -0,0 +1,71 @@
+#define PY_SSIZE_T_CLEAN
+#include <Python.h>
+
+int
+main(int argc, char *argv[])
+{
+    PyObject *pName, *pModule, *pFunc;
+    PyObject *pArgs, *pValue;
+    int i;
+
+    if (argc < 3) {
+        fprintf(stderr,"Usage: call pythonfile funcname [args]\n");
+        return 1;
+    }
+
+    Py_Initialize();
+    pName = PyUnicode_DecodeFSDefault(argv[1]);
+    /* Error checking of pName left out */
+    pModule = PyImport_Import(pName);
+
+    Py_DECREF(pName);
+
+    if (pModule != NULL) {
+        pFunc = PyObject_GetAttrString(pModule, argv[2]);
+        /* pFunc is a new reference */
+
+        if (pFunc && PyCallable_Check(pFunc)) {
+            pArgs = PyTuple_New(argc - 3);
+            for (i = 0; i < argc - 3; ++i) {
+                pValue = PyLong_FromLong(atoi(argv[i + 3]));
+                if (!pValue) {
+                    Py_DECREF(pArgs);
+                    Py_DECREF(pModule);
+                    fprintf(stderr, "Cannot convert argument\n");
+                    return 1;
+                }
+                /* pValue reference stolen here: */
+                PyTuple_SetItem(pArgs, i, pValue);
+            }
+            pValue = PyObject_CallObject(pFunc, pArgs);
+            Py_DECREF(pArgs);
+            if (pValue != NULL) {
+                printf("Result of call: %ld\n", PyLong_AsLong(pValue));
+                Py_DECREF(pValue);
+            }
+            else {
+                Py_DECREF(pFunc);
+                Py_DECREF(pModule);
+                PyErr_Print();
+                fprintf(stderr,"Call failed\n");
+                return 1;
+            }
+        }
+        else {
+            if (PyErr_Occurred())
+                PyErr_Print();
+            fprintf(stderr, "Cannot find function \"%s\"\n", argv[2]);
+        }
+        Py_XDECREF(pFunc);
+        Py_DECREF(pModule);
+    }
+    else {
+        PyErr_Print();
+        fprintf(stderr, "Failed to load \"%s\"\n", argv[1]);
+        return 1;
+    }
+    if (Py_FinalizeEx() < 0) {
+        return 120;
+    }
+    return 0;
+}
diff --git a/embed/mult.py b/embed/mult.py
new file mode 100644
index 0000000..97ca688
--- /dev/null
+++ b/embed/mult.py
@@ -0,0 +1,7 @@
+def multiply(a,b):
+    print("Will compute", a, "times", b)
+    c = 0
+    for i in range(0, a):
+        c = c + b
+    return c
+