diff --git a/.gitignore b/.gitignore
index a57a1904..0932e2b8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,6 @@ robotDescription.cpp
 log_files/
 
 *Example
+*.so
+*.egg-info
+__pycache__
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f4c4ac27..9f68bd3d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -50,6 +50,14 @@ else ()
     message(STATUS "Building with Magnum disabled")
 endif ()
 
+option(WITH_PYBIND OFF)
+if (WITH_PYBIND)
+    find_package(pybind11 REQUIRED)
+else ()
+    message(STATUS "Building with pybind disabled")
+endif ()
+
+
 if (${HAVE_SYMENGINE_LLVM})
     message(STATUS "Testing SymEngine LLVM & SBML support - found")
 else ()
@@ -199,6 +207,18 @@ else()
     message(WARNING "No robotDescription.cpp found. No executable will be created.")
 endif()
 
+if (WITH_PYBIND)
+    pybind11_add_module(py_dismech
+            src/app.cpp
+    )
+    target_link_libraries(py_dismech PRIVATE
+            common_sources
+    )
+
+    set_target_properties(py_dismech PROPERTIES
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/py_dismech"
+)
+endif()
 
 
 option(CREATE_EXAMPLES OFF)
diff --git a/README.md b/README.md
index 102359ab..1b218226 100755
--- a/README.md
+++ b/README.md
@@ -41,8 +41,10 @@ Based on the [Discrete Elastic Rods](https://www.cs.columbia.edu/cg/pdfs/143-rod
 ### TODO
 If you'd like DisMech to support a new feature, feel free create an issue and we'll add it to the list here.
 
+#### Ongoing
+- [ ] Expand Python bindings. PR [#13](https://github.com/StructuresComp/dismech-rods/pull/13)
+
 #### High priority
-- [ ] Migrate API to Python via pybind11
 - [ ] Add extension control via natural length manipulation.
 - [ ] Add shell functionality.
 - [ ] Improve robustness of friction.
@@ -177,10 +179,31 @@ Below is an example rendering.
     make -j$(nproc)
     ```
   - For those wishing to customize the rendering settings, users can take a look and adjust the source code in `magnumDERSimulationEnvironment.cpp` accordingly.
+- [Pybind11](https://github.com/pybind/pybind11)
+  - DisMech also offers Python bindings via pybind11. Users can either install by source or via `pip install pybind11`.
+  Afterward, the bindings can be built as follows:
+    ```bash
+    mkdir build && cd build
+    cmake -DWITH_PYBIND=on ..
+    make -j$(nproc)
+    cd ..
+    pip install -e .  # this will install py_dismech
+    ```
+  - Note that when building with older versions of `MKL` with Python, it may be possible to run into the following [static linkage bug](https://bugs.launchpad.net/ubuntu/+source/intel-mkl/+bug/1947626).
+  Users can either upgrade their `MKL` version or specify their `LD_PRELOAD` explicitly to include the following .so files.
+    ```bash
+    export LD_PRELOAD=$MKL_LIB/libmkl_def.so.2:\
+                      $MKL_LIB/libmkl_avx2.so.2:\
+                      $MKL_LIB/libmkl_core.so.2:\
+                      $MKL_LIB/libmkl_intel_lp64.so.2:\
+                      $MKL_LIB/libmkl_intel_thread.so.2:\
+                      /usr/lib/x86_64-linux-gnu/libiomp5.so
+     ```
+  - All available bindings can be seen in `app.cpp`. Users can expect heavy development for expanding the python bindings as time goes on.
 
 ***
 
-### Running Examples
+### Running Examples in C++
 DisMech is setup so that simulation environments can be instantiated using a single cpp file called `robotDescription.cpp`.
 
 Several example of working DisMech simulations can be seen in the `examples/` directory.
@@ -214,6 +237,12 @@ in `/examples`, their Jacobian matrices are small enough that any amount of para
 simulation. Therefore, it is recommended to set `OMP_NUM_THREADS=1` (`dismech.sh` does this as automatically) and see if parallelization is worth it
 for larger systems through profiling.
 
+### Running Examples in Python
+
+Python versions of certain examples can be found in `py_examples/`. To run, simply run the provided scripts:
+```bash
+OMP_NUM_THREADS=1 python py_examples/spider_case/spiderExample.py
+```
 ***
 
 ### Creating Custom Simulation Environments
diff --git a/py_dismech/__init__.py b/py_dismech/__init__.py
new file mode 100644
index 00000000..d9db50a6
--- /dev/null
+++ b/py_dismech/__init__.py
@@ -0,0 +1 @@
+from .py_dismech import *
diff --git a/py_examples/helix_case/helixExample.py b/py_examples/helix_case/helixExample.py
new file mode 100644
index 00000000..7eac26cc
--- /dev/null
+++ b/py_examples/helix_case/helixExample.py
@@ -0,0 +1,56 @@
+"""
+    Helix Under Gravity Example
+
+    This is a pybind version of helixExample.cpp
+
+"""
+import numpy as np
+import sys
+import py_dismech
+from pathlib import Path
+from functools import partial
+
+sim_manager = py_dismech.SimulationManager()
+
+soft_robots = sim_manager.soft_robots
+sim_params = sim_manager.sim_params
+render_params = sim_manager.render_params
+create_joint = sim_manager.soft_robots.createJoint
+add_to_joint = sim_manager.soft_robots.addToJoint
+add_force = sim_manager.forces.addForce
+
+############################
+
+sim_params.dt = 5e-3
+sim_params.sim_time = 10
+sim_params.dtol = 1e-3
+sim_params.integrator = py_dismech.IMPLICIT_MIDPOINT
+
+render_params.render_scale = 5.0
+render_params.render_per = 5
+
+# Read vertices describing helical shape from a file
+vertices = np.loadtxt(Path(__file__).parents[2] / 'examples/helix_case/helix_configuration.txt')
+
+# Create the helix limb with custom configuration
+radius = 5e-3
+young_mod = 1e7
+density = 1273.52
+poisson = 0.5
+
+add_limb = partial(sim_manager.soft_robots.addLimb, rho=density, rod_radius=radius,
+                   youngs_modulus=young_mod, poisson_ratio=poisson)
+
+# Add the helical structure as a sequential series of vertices
+add_limb(vertices)
+
+# Fix the top end of the helix (locking the first node)
+soft_robots.lockEdge(0, 0)
+
+# Add gravity
+gravity_force = py_dismech.GravityForce(soft_robots, np.array([0.0, 0.0, -9.8]))
+add_force(gravity_force)
+
+# Initialize and run the simulation
+sim_manager.initialize(sys.argv)
+sim_manager.run()
diff --git a/py_examples/spider_case/spiderExample.py b/py_examples/spider_case/spiderExample.py
new file mode 100644
index 00000000..1ff41e09
--- /dev/null
+++ b/py_examples/spider_case/spiderExample.py
@@ -0,0 +1,88 @@
+"""
+    Spider on Incline Example
+
+    This is a pybind version of spiderExample.cpp
+"""
+
+import numpy as np
+import sys
+import py_dismech
+from functools import partial
+
+sim_manager = py_dismech.SimulationManager()
+
+soft_robots = sim_manager.soft_robots
+sim_params = sim_manager.sim_params
+render_params = sim_manager.render_params
+create_joint = sim_manager.soft_robots.createJoint
+add_to_joint = sim_manager.soft_robots.addToJoint
+add_force = sim_manager.forces.addForce
+
+SIM_FAST = True
+
+############################
+
+sim_params.sim_time = 2
+sim_params.ftol = 1e-3
+
+render_params.render_scale = 5.0
+render_params.show_mat_frames = True
+
+if SIM_FAST:
+    sim_params.dt = 2.5e-3
+    sim_params.max_iter.num_iters = 15
+    sim_params.max_iter.terminate_at_max = False
+    delta = 5e-3
+    nu = 1e-2
+else:
+    sim_params.dt = 1e-3
+    sim_params.adaptive_time_stepping = 7
+    delta = 5e-4
+    nu = 5e-3
+
+
+n = 25
+radius = 5e-3
+young_mod = 3e6
+density = 1180
+poisson = 0.5
+mu = 0.4
+add_limb = partial(sim_manager.soft_robots.addLimb, num_nodes=n, rho=density, rod_radius=radius,
+                   youngs_modulus=young_mod, poisson_ratio=poisson, mu=mu)
+
+add_limb(np.array([ 0.00,  0.00,  0.20]),    np.array([ 0.00,  0.00,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([ 0.10,  0.00,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([ 0.00,  0.10,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([ 0.00, -0.10,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([-0.10,  0.00,  0.10]))
+add_limb(np.array([ 0.10,  0.00,  0.10]),    np.array([ 0.10,  0.00,  0.00]))
+add_limb(np.array([ 0.00,  0.10,  0.10]),    np.array([ 0.00,  0.10,  0.00]))
+add_limb(np.array([ 0.00, -0.10,  0.10]),    np.array([ 0.00, -0.10,  0.00]))
+add_limb(np.array([-0.10,  0.00,  0.10]),    np.array([-0.10,  0.00,  0.00]))
+
+create_joint(0, -1)
+add_to_joint(0, 1, 0)
+add_to_joint(0, 2, 0)
+add_to_joint(0, 3, 0)
+add_to_joint(0, 4, 0)
+create_joint(1, -1)
+add_to_joint(1, 5, 0)
+create_joint(2, -1)
+add_to_joint(2, 6, 0)
+create_joint(3, -1)
+add_to_joint(3, 7, 0)
+create_joint(4, -1)
+add_to_joint(4, 8, 0)
+
+# Add gravity with a slight x-axis perturbation
+gravity_force = py_dismech.GravityForce(soft_robots, np.array([1.0, 0.0, -9.8]))
+add_force(gravity_force)
+
+# Add floor contact
+floor_z = -0.10
+floor_contact_force = py_dismech.FloorContactForce(soft_robots, delta, nu, floor_z)
+add_force(floor_contact_force)
+
+# Initialize and run the simulation
+sim_manager.initialize(sys.argv)
+sim_manager.run()
diff --git a/setup.py b/setup.py
new file mode 100644
index 00000000..a1dce66d
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,14 @@
+from setuptools import setup, Extension
+import glob
+import os
+
+# Use glob to find the actual .so file with cp310-... in its name
+so_file = glob.glob(os.path.join('py_dismech', 'py_dismech*.so'))[0]  # Adjust path as needed
+
+setup(
+    name='py_dismech',
+    version='0.1',
+    packages=['py_dismech'],
+    package_dir={'py_dismech': 'py_dismech'},
+    package_data={'py_dismech': [os.path.basename(so_file)]},  # Include the .so file
+)
diff --git a/src/app.cpp b/src/app.cpp
new file mode 100644
index 00000000..b5e2f31d
--- /dev/null
+++ b/src/app.cpp
@@ -0,0 +1,185 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/functional.h>
+#include <pybind11/numpy.h>
+#include <pybind11/eigen.h>
+#include <pybind11/stl.h>
+
+#include "rod_mechanics/external_forces/dampingForce.h"
+#include "rod_mechanics/external_forces/gravityForce.h"
+#include "rod_mechanics/external_forces/floorContactForce.h"
+#include "rod_mechanics/external_forces/uniformConstantForce.h"
+#include "rod_mechanics/external_forces/contactForce.h"
+
+#include "globalDefinitions.h"
+#include "world.h"
+#include "logging/worldLogger.h"
+#include "simulation_environments/headlessDERSimulationEnvironment.h"
+#include "simulation_environments/openglDERSimulationEnvironment.h"
+#ifdef WITH_MAGNUM
+#include "simulation_environments/magnumDERSimulationEnvironment.h"
+#endif
+
+
+namespace py = pybind11;
+
+
+// Hack: main creates the output file for logging
+ofstream logging_output_file;
+
+double openglDERSimulationEnvironment::render_scale = 1.0;
+bool openglDERSimulationEnvironment::show_mat_frames = false;
+
+// Forward cmdline args from python
+int argc;
+std::vector<char*> argv;
+
+
+// Wrapper for running the simulation from Python
+class simulationManager {
+public:
+    std::shared_ptr<softRobots> soft_robots;
+    std::shared_ptr<forceContainer> forces;
+    std::shared_ptr<world> my_world;
+    simParams sim_params;
+    renderParams render_params;
+    std::shared_ptr<worldLogger> logger;
+    std::unique_ptr<derSimulationEnvironment> env;
+
+    simulationManager() {
+        soft_robots = std::make_shared<softRobots>();
+        forces = std::make_shared<forceContainer>();
+        logger = nullptr;
+    }
+
+    void initialize(int argc, char* argv[]) {
+        soft_robots->setup();
+        my_world = std::make_shared<world>(soft_robots, forces, sim_params);
+
+        switch (render_params.renderer) {
+            case HEADLESS:
+                env = std::make_unique<headlessDERSimulationEnvironment>(my_world, render_params, logger);
+                break;
+            case OPENGL:
+                env = std::make_unique<openglDERSimulationEnvironment>(my_world, render_params, logger, argc, argv);
+                break;
+#ifdef WITH_MAGNUM
+            case MAGNUM:
+                env = std::make_unique<Magnum::magnumDERSimulationEnvironment>(my_world, render_params, logger, argc, argv);
+                break;
+#endif
+            default:
+                throw std::runtime_error("Unknown renderer type provided.");
+        }
+    }
+
+    void run() {
+        if (env) {
+            env->runSimulation();
+        }
+    }
+};
+
+
+PYBIND11_MODULE(py_dismech, m) {
+    // Wrapper class for accessing dismech data and logic via pybind
+    py::class_<simulationManager>(m, "SimulationManager")
+        .def(py::init<>())
+        .def("initialize", [](simulationManager& self, std::vector<std::string> args) {
+            // Convert std::vector<std::string> to argc and argv
+            argc = args.size();
+            for (int i = 0; i < argc; ++i) {
+                argv.push_back(const_cast<char*>(args[i].c_str()));
+            }
+            argv.push_back(nullptr);  // Null-terminate the argv array
+            // Call the original function
+            self.initialize(argc, argv.data());
+        })
+        .def("run", &simulationManager::run)
+        .def_readwrite("soft_robots", &simulationManager::soft_robots)
+        .def_readwrite("forces", &simulationManager::forces)
+        .def_readwrite("sim_params", &simulationManager::sim_params)
+        .def_readwrite("render_params", &simulationManager::render_params)
+        .def_readwrite("logger", &simulationManager::logger);
+
+    py::class_<softRobots, std::shared_ptr<softRobots>>(m, "SoftRobots")
+        .def(py::init<>())
+        .def("addLimb",
+            py::overload_cast<const Eigen::Vector3d&, const Eigen::Vector3d&,
+                              int, double, double, double, double, double>(&softRobots::addLimb),
+            py::arg("start"), py::arg("end"), py::arg("num_nodes"),
+            py::arg("rho"), py::arg("rod_radius"), py::arg("youngs_modulus"),
+            py::arg("poisson_ratio"), py::arg("mu") = 0.0)
+        .def("addLimb",
+            py::overload_cast<const std::vector<Eigen::Vector3d>&, double, double,
+                              double, double, double>(&softRobots::addLimb),
+            py::arg("nodes"), py::arg("rho"), py::arg("rod_radius"),
+            py::arg("youngs_modulus"), py::arg("poisson_ratio"), py::arg("mu") = 0.0)
+        .def("createJoint", &softRobots::createJoint, py::arg("limb_idx"), py::arg("node_idx"))
+        .def("addToJoint", &softRobots::addToJoint, py::arg("joint_idx"), py::arg("limb_idx"), py::arg("node_idx"))
+        .def("lockEdge", &softRobots::lockEdge, py::arg("limb_idx"), py::arg("edge_idx"))
+        .def("applyInitialVelocities", &softRobots::applyInitialVelocities, py::arg("limb_idx"), py::arg("velocities"))
+        .def("setup", &softRobots::setup)
+        .def("addController", &softRobots::addController, py::arg("controller"))
+        .def_readwrite("limbs", &softRobots::limbs)
+        .def_readwrite("joints", &softRobots::joints)
+        .def_readwrite("controllers", &softRobots::controllers);
+
+    // ====================================== Enum Definitions ============================================
+    py::enum_<integratorMethod>(m, "IntegratorMethod")
+        .value("FORWARD_EULER", FORWARD_EULER)
+        .value("VERLET_POSITION", VERLET_POSITION)
+        .value("BACKWARD_EULER", BACKWARD_EULER)
+        .value("IMPLICIT_MIDPOINT", IMPLICIT_MIDPOINT)
+        .export_values();
+
+    py::enum_<renderEngine>(m, "RenderEngine")
+        .value("HEADLESS", HEADLESS)
+        .value("OPENGL", OPENGL)
+        .value("MAGNUM", MAGNUM)
+        .export_values();
+
+    // ===================================== Struct Definitions ============================================
+    py::class_<simParams::maxIterations>(m, "MaxIterations")
+        .def(py::init<>())
+        .def_readwrite("num_iters", &simParams::maxIterations::num_iters)
+        .def_readwrite("terminate_at_max", &simParams::maxIterations::terminate_at_max);
+
+    py::class_<simParams>(m, "SimParams")
+        .def(py::init<>())
+        .def_readwrite("sim_time", &simParams::sim_time)
+        .def_readwrite("dt", &simParams::dt)
+        .def_readwrite("integrator", &simParams::integrator)
+        .def_readwrite("dtol", &simParams::dtol)
+        .def_readwrite("ftol", &simParams::ftol)
+        .def_readwrite("max_iter", &simParams::max_iter)
+        .def_readwrite("line_search", &simParams::line_search)
+        .def_readwrite("adaptive_time_stepping", &simParams::adaptive_time_stepping)
+        .def_readwrite("enable_2d_sim", &simParams::enable_2d_sim);
+
+    py::class_<renderParams>(m, "RenderParams")
+        .def(py::init<>())
+        .def_readwrite("renderer", &renderParams::renderer)
+        .def_readwrite("render_scale", &renderParams::render_scale)
+        .def_readwrite("cmd_line_per", &renderParams::cmd_line_per)
+        .def_readwrite("render_per", &renderParams::render_per)
+        .def_readwrite("render_record_path", &renderParams::render_record_path)
+        .def_readwrite("show_mat_frames", &renderParams::show_mat_frames)
+        .def_readwrite("debug_verbosity", &renderParams::debug_verbosity);
+
+    // ====================================== Force Definitions ============================================
+    py::class_<forceContainer, std::shared_ptr<forceContainer>>(m, "ForceContainer")
+        .def(py::init<>())
+        .def(py::init<const std::vector<std::shared_ptr<baseForce>>&>(), py::arg("m_forces"))
+        .def("addForce", &forceContainer::addForce, py::arg("force"));
+
+    py::class_<baseForce, std::shared_ptr<baseForce>>(m, "BaseForce");
+
+    py::class_<gravityForce, std::shared_ptr<gravityForce>, baseForce>(m, "GravityForce")
+        .def(py::init<const std::shared_ptr<softRobots>&, const Eigen::Vector3d &>(),
+            py::arg("soft_robots"), py::arg("g_vector"));
+
+    py::class_<floorContactForce, std::shared_ptr<floorContactForce>, baseForce>(m, "FloorContactForce")
+        .def(py::init<const std::shared_ptr<softRobots>&, double, double, double, double>(),
+            py::arg("soft_robots"), py::arg("floor_delta"), py::arg("floor_slipTol"),
+            py::arg("floor_z"), py::arg("floor_mu") = 0.0);
+}