diff --git a/src/extern/include/gridfire/extern/gridfire_context.h b/src/extern/include/gridfire/extern/gridfire_context.h
new file mode 100644
index 00000000..5c432820
--- /dev/null
+++ b/src/extern/include/gridfire/extern/gridfire_context.h
@@ -0,0 +1,40 @@
+#ifndef GF_GRIDFIRE_CONTEXT_H
+#define GF_GRIDFIRE_CONTEXT_H
+
+#include "gridfire/gridfire.h"
+#include "fourdst/atomic/atomicSpecies.h"
+
+#include <memory>
+#include <vector>
+
+struct GridFireContext {
+    std::unique_ptr<gridfire::policy::NetworkPolicy> policy;
+    gridfire::engine::DynamicEngine* engine;
+    std::unique_ptr<gridfire::solver::DynamicNetworkSolverStrategy> solver;
+
+    std::vector<fourdst::atomic::Species> speciesList;
+    fourdst::composition::Composition working_comp;
+
+    void init_species_map(const std::vector<std::string>& species_names);
+    void init_engine_from_policy(const std::string& policy_name, const double *abundances, size_t num_species);
+    void init_solver_from_engine(const std::string& solver_name);
+
+    void init_composition_from_abundance_vector(const double* abundances, size_t num_species);
+
+    int evolve(
+        const double* Y_in,
+        size_t num_species,
+        double T,
+        double rho,
+        double tMax,
+        double dt0,
+        double* Y_out,
+        double& energy_out,
+        double& dEps_dT,
+        double& dEps_dRho, double& mass_lost
+    );
+
+    std::string last_error;
+};
+
+#endif
\ No newline at end of file
diff --git a/src/extern/include/gridfire/extern/gridfire_extern.h b/src/extern/include/gridfire/extern/gridfire_extern.h
new file mode 100644
index 00000000..5c535fc8
--- /dev/null
+++ b/src/extern/include/gridfire/extern/gridfire_extern.h
@@ -0,0 +1,91 @@
+#ifndef GF_GRIDFIRE_EXTERN_H
+#define GF_GRIDFIRE_EXTERN_H
+
+#include <stddef.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+    enum FDSSE_ERROR_CODES {
+        FDSSE_NON_4DSTAR_ERROR = -102,
+        FDSSE_UNKNOWN_ERROR = -101,
+        FDSSE_SUCCESS = 1,
+        FDSSE_UNKNOWN_SYMBOL_ERROR = 100,
+        FDSSE_SPECIES_ERROR = 101,
+        FDSSE_INVALID_COMPOSITION_ERROR = 102,
+        FDSSE_COMPOSITION_ERROR = 103
+    };
+
+    enum GF_ERROR_CODES {
+        GF_NON_GRIDFIRE_ERROR = -2,
+        GF_UNKNOWN_ERROR = -1,
+        GF_SUCCESS = 0,
+
+        GF_INVALID_QSE_SOLUTION_ERROR = 5,
+        GF_FAILED_TO_PARTITION_ENGINE_ERROR = 6,
+        GF_NETWORK_RESIZED_ERROR = 7,
+        GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR = 8,
+        GF_BAD_COLLECTION_ERROR = 9,
+        GF_BAD_RHS_ENGINE_ERROR = 10,
+        GF_STALE_JACOBIAN_ERROR = 11,
+        GF_UNINITIALIZED_JACOBIAN_ERROR = 12,
+        GF_UNKNOWN_JACOBIAN_ERROR = 13,
+        GF_JACOBIAN_ERROR = 14,
+        GF_ENGINE_ERROR = 15,
+
+        GF_MISSING_BASE_REACTION_ERROR = 16,
+        GF_MISSING_SEED_SPECIES_ERROR = 17,
+        GF_MISSING_KEY_REACTION_ERROR = 18,
+        GF_POLICY_ERROR = 19,
+
+        GF_REACTION_PARSING_ERROR = 20,
+        GF_REACTION_ERROR = 21,
+
+        GF_SINGULAR_JACOBIAN_ERROR = 22,
+        GF_ILL_CONDITIONED_JACOBIAN_ERROR = 23,
+        GF_CVODE_SOLVER_FAILURE_ERROR = 24,
+        GF_KINSOL_SOLVER_FAILURE_ERROR = 25,
+        GF_SUNDIALS_ERROR = 26,
+        GF_SOLVER_ERROR = 27,
+
+        GF_HASHING_ERROR = 28,
+        GF_UTILITY_ERROR = 29,
+
+        GF_DEBUG_ERROR = 30,
+
+        GF_GRIDFIRE_ERROR = 31,
+    };
+
+    char* gf_get_last_error_message(void* ptr);
+
+    char* gf_error_code_to_string(int error_code);
+
+    void* gf_init();
+
+    int gf_free(void* ctx);
+
+    int gf_register_species(void* ptr, const int num_species, const char** species_names);
+
+    int gf_construct_engine_from_policy(void* ptr, const char* policy_name, const double *abundances, size_t num_species);
+
+    int gf_construct_solver_from_engine(void* ptr, const char* solver_name);
+
+    int gf_evolve(
+        void* ptr,
+        const double* Y_in,
+        size_t num_species,
+        double T,
+        double rho,
+        double tMax,
+        double dt0,
+        double* Y_out,
+        double* energy_out,
+        double* dEps_dT,
+        double* dEps_dRho, double* mass_lost
+    );
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/src/extern/lib/gridfire_context.cpp b/src/extern/lib/gridfire_context.cpp
new file mode 100644
index 00000000..4df5c68f
--- /dev/null
+++ b/src/extern/lib/gridfire_context.cpp
@@ -0,0 +1,156 @@
+#include "gridfire/extern/gridfire_context.h"
+
+#include "fourdst/atomic/species.h"
+
+#include "fourdst/composition/exceptions/exceptions_composition.h"
+
+void GridFireContext::init_species_map(const std::vector<std::string> &species_names) {
+    for (const auto& name: species_names) {
+        working_comp.registerSymbol(name);
+    }
+    this->speciesList.clear();
+    this->speciesList.reserve(species_names.size());
+
+    auto resolve_species_name = [](const std::string& name) -> fourdst::atomic::Species {
+        if (fourdst::atomic::species.contains(name)) {
+            return fourdst::atomic::species.at(name);
+        }
+        throw fourdst::composition::exceptions::UnknownSymbolError("Species " + name + " is not recognized in the atomic species database.");
+    };
+
+    for (const auto& name: species_names) {
+        this->speciesList.push_back(resolve_species_name(name));
+    }
+
+}
+
+void GridFireContext::init_engine_from_policy(const std::string &policy_name, const double *abundances, const size_t num_species) {
+    init_composition_from_abundance_vector(abundances, num_species);
+
+    enum class EnginePolicy {
+        MAIN_SEQUENCE_POLICY
+    };
+
+    static const std::unordered_map<std::string, EnginePolicy> engine_map = {
+        {"MAIN_SEQUENCE_POLICY", EnginePolicy::MAIN_SEQUENCE_POLICY}
+    };
+
+    if (!engine_map.contains(policy_name)) {
+        throw gridfire::exceptions::PolicyError(
+            std::format(
+                "Engine Policy {} is not recognized. Valid policies are: {}",
+                policy_name,
+                gridfire::utils::iterable_to_delimited_string(engine_map, ", ", [](const auto& pair){ return pair.first; })
+            )
+        );
+    }
+
+    switch (engine_map.at(policy_name)) {
+        case EnginePolicy::MAIN_SEQUENCE_POLICY: {
+            this->policy = std::make_unique<gridfire::policy::MainSequencePolicy>(this->working_comp);
+            this->engine = &policy->construct();
+            break;
+        }
+        default:
+            throw gridfire::exceptions::PolicyError(
+                "Unhandled engine policy in GridFireContext::init_engine_from_policy"
+            );
+    }
+
+
+}
+
+void GridFireContext::init_solver_from_engine(const std::string &solver_name) {
+    enum class SolverType {
+        CVODE
+    };
+
+    static const std::unordered_map<std::string, SolverType> solver_map = {
+        {"CVODE", SolverType::CVODE}
+    };
+
+    if (!solver_map.contains(solver_name)) {
+        throw gridfire::exceptions::SolverError(
+            std::format(
+                "Solver {} is not recognized. Valid solvers are: {}",
+                solver_name,
+                gridfire::utils::iterable_to_delimited_string(solver_map, ", ", [](const auto& pair){ return pair.first; })
+            )
+        );
+    }
+
+    switch (solver_map.at(solver_name)) {
+        case SolverType::CVODE: {
+            this->solver = std::make_unique<gridfire::solver::CVODESolverStrategy>(*this->engine);
+            break;
+        }
+        default:
+            throw gridfire::exceptions::SolverError(
+                "Unhandled solver type in GridFireContext::init_solver_from_engine"
+            );
+    }
+
+}
+
+void GridFireContext::init_composition_from_abundance_vector(const double *abundances, size_t num_species) {
+    if (num_species == 0) {
+        throw fourdst::composition::exceptions::InvalidCompositionError("Cannot initialize composition with zero species.");
+    }
+    if (num_species != working_comp.size()) {
+        throw fourdst::composition::exceptions::InvalidCompositionError(
+            std::format(
+                "Number of species provided ({}) does not match the registered species count ({}).",
+                num_species,
+                working_comp.size()
+            )
+        );
+    }
+
+    for (size_t i = 0; i < num_species; i++) {
+        this->working_comp.setMolarAbundance(this->speciesList[i], abundances[i]);
+    }
+}
+
+int GridFireContext::evolve(
+    const double* Y_in,
+    const size_t num_species,
+    const double T,
+    const double rho,
+    const double tMax,
+    const double dt0,
+    double* Y_out,
+    double& energy_out,
+    double& dEps_dT,
+    double& dEps_dRho, double& mass_lost
+) {
+    init_composition_from_abundance_vector(Y_in, num_species);
+
+    gridfire::NetIn netIn;
+    netIn.temperature = T;
+    netIn.density = rho;
+    netIn.dt0 = dt0;
+    netIn.tMax = tMax;
+    netIn.composition = this->working_comp;
+
+    const gridfire::NetOut result = this->solver->evaluate(netIn);
+
+    energy_out = result.energy;
+    dEps_dT  = result.dEps_dT;
+    dEps_dRho = result.dEps_dRho;
+
+    std::set<fourdst::atomic::Species> seen_species;
+    for (size_t i = 0; i < num_species; i++) {
+        fourdst::atomic::Species species = this->speciesList[i];
+        Y_out[i] = result.composition.getMolarAbundance(species);
+        seen_species.insert(species);
+    }
+
+    mass_lost = 0.0;
+    for (const auto& species : result.composition.getRegisteredSpecies()) {
+        if (!seen_species.contains(species)) {
+            mass_lost += species.mass() * result.composition.getMolarAbundance(species);
+        }
+    }
+
+    return 0;
+}
\ No newline at end of file
diff --git a/src/extern/lib/gridfire_extern.cpp b/src/extern/lib/gridfire_extern.cpp
new file mode 100644
index 00000000..b51a7edb
--- /dev/null
+++ b/src/extern/lib/gridfire_extern.cpp
@@ -0,0 +1,294 @@
+#include "gridfire/gridfire.h"
+#include "fourdst/composition/exceptions/exceptions_composition.h"
+#include "gridfire/extern/gridfire_context.h"
+#include "gridfire/extern/gridfire_extern.h"
+
+extern "C" {
+    void* gf_init() {
+        return new GridFireContext();
+    }
+
+    int gf_free(void* ctx) {
+        delete static_cast<GridFireContext*>(ctx);
+        return 0;
+    }
+
+    int gf_register_species(void* ptr, const int num_species, const char** species_names) {
+        auto* ctx = static_cast<GridFireContext*>(ptr);
+        try {
+            std::vector<std::string> names;
+            for(int i=0; i<num_species; ++i) {
+                names.emplace_back(species_names[i]);
+            }
+            ctx->init_species_map(names);
+            return FDSSE_SUCCESS;
+        } catch (const fourdst::composition::exceptions::UnknownSymbolError& e) {
+            ctx->last_error = e.what();
+            return FDSSE_UNKNOWN_SYMBOL_ERROR;
+        } catch (const fourdst::composition::exceptions::SpeciesError& e) {
+            ctx->last_error = e.what();
+            return FDSSE_SPECIES_ERROR;
+        } catch (const std::exception& e) {
+            ctx->last_error = e.what();
+            return FDSSE_NON_4DSTAR_ERROR;
+        } catch (...) {
+            ctx->last_error = "Unknown error occurred during species registration.";
+            return FDSSE_UNKNOWN_ERROR;
+        }
+    }
+
+    int gf_construct_engine_from_policy(
+        void* ptr,
+        const char* policy_name,
+        const double *abundances,
+        const size_t num_species
+    ) {
+        auto* ctx = static_cast<GridFireContext*>(ptr);
+        try {
+            ctx->init_engine_from_policy(std::string(policy_name), abundances, num_species);
+            return GF_SUCCESS;
+        } catch (const gridfire::exceptions::MissingBaseReactionError& e) {
+            ctx->last_error = e.what();
+            return GF_MISSING_BASE_REACTION_ERROR;
+        } catch (const gridfire::exceptions::MissingSeedSpeciesError& e) {
+            ctx->last_error = e.what();
+            return GF_MISSING_SEED_SPECIES_ERROR;
+        }  catch (const gridfire::exceptions::MissingKeyReactionError& e) {
+            ctx->last_error = e.what();
+            return GF_MISSING_KEY_REACTION_ERROR;
+        } catch (const gridfire::exceptions::PolicyError& e) {
+            ctx->last_error = e.what();
+            return GF_POLICY_ERROR;
+        } catch (std::exception& e) {
+            ctx->last_error = e.what();
+            return GF_NON_GRIDFIRE_ERROR;
+        } catch (...) {
+            ctx->last_error = "Unknown error occurred during engine construction.";
+            return GF_UNKNOWN_ERROR;
+        }
+    }
+
+    int gf_construct_solver_from_engine(
+        void* ptr,
+        const char* solver_name
+    ) {
+        auto* ctx = static_cast<GridFireContext*>(ptr);
+        try {
+            ctx->init_solver_from_engine(std::string(solver_name));
+            return GF_SUCCESS;
+        } catch (std::exception& e) {
+            ctx->last_error = e.what();
+            return GF_NON_GRIDFIRE_ERROR;
+        } catch (...) {
+            ctx->last_error = "Unknown error occurred during solver construction.";
+            return GF_UNKNOWN_ERROR;
+        }
+    }
+
+    int gf_evolve(
+        void* ptr,
+        const double* Y,
+        const size_t num_species,
+        const double T,
+        const double rho,
+        const double tMax,
+        const double dt0,
+        double* Y_out,
+        double* energy_out,
+        double* dEps_dT,
+        double* dEps_dRho, double* mass_lost
+    ) {
+        auto* ctx = static_cast<GridFireContext*>(ptr);
+        try {
+            const int result = ctx->evolve(
+                Y,
+                num_species,
+                T,
+                rho,
+                tMax,
+                dt0,
+                Y_out,
+                *energy_out,
+                *dEps_dT,
+                *dEps_dRho, *mass_lost
+            );
+            if (result != 0) {
+                return result;
+            }
+            return GF_SUCCESS;
+        } catch (fourdst::composition::exceptions::UnknownSymbolError& e) {
+            ctx->last_error = e.what();
+            return FDSSE_UNKNOWN_SYMBOL_ERROR;
+        } catch (const fourdst::composition::exceptions::SpeciesError& e) {
+            ctx->last_error = e.what();
+            return FDSSE_SPECIES_ERROR;
+        } catch (const fourdst::composition::exceptions::InvalidCompositionError& e) {
+            ctx->last_error = e.what();
+            return FDSSE_INVALID_COMPOSITION_ERROR;
+        } catch (const fourdst::composition::exceptions::CompositionError& e) {
+            ctx->last_error = e.what();
+            return FDSSE_COMPOSITION_ERROR;
+        } catch (const gridfire::exceptions::InvalidQSESolutionError& e) {
+            ctx->last_error = e.what();
+            return GF_INVALID_QSE_SOLUTION_ERROR;
+        } catch (const gridfire::exceptions::FailedToPartitionEngineError& e) {
+            ctx->last_error = e.what();
+            return GF_FAILED_TO_PARTITION_ENGINE_ERROR;
+        } catch (const gridfire::exceptions::NetworkResizedError& e) {
+            ctx->last_error = e.what();
+            return GF_NETWORK_RESIZED_ERROR;
+        } catch (const gridfire::exceptions::UnableToSetNetworkReactionsError& e) {
+            ctx->last_error = e.what();
+            return GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR;
+        } catch (const gridfire::exceptions::BadCollectionError& e) {
+            ctx->last_error = e.what();
+            return GF_BAD_COLLECTION_ERROR;
+        } catch (const gridfire::exceptions::BadRHSEngineError& e) {
+            ctx->last_error = e.what();
+            return GF_BAD_RHS_ENGINE_ERROR;
+        } catch (const gridfire::exceptions::StaleJacobianError& e) {
+            ctx->last_error = e.what();
+            return GF_STALE_JACOBIAN_ERROR;
+        } catch (const gridfire::exceptions::UninitializedJacobianError& e) {
+            ctx->last_error = e.what();
+            return GF_UNINITIALIZED_JACOBIAN_ERROR;
+        } catch (const gridfire::exceptions::UnknownJacobianError& e) {
+            ctx->last_error = e.what();
+            return GF_UNKNOWN_JACOBIAN_ERROR;
+        } catch (const gridfire::exceptions::JacobianError& e) {
+            ctx->last_error = e.what();
+            return GF_JACOBIAN_ERROR;
+        } catch (const gridfire::exceptions::EngineError& e) {
+            ctx->last_error = e.what();
+            return GF_ENGINE_ERROR;
+        } catch (const gridfire::exceptions::ReactionParsingError& e) {
+            ctx->last_error = e.what();
+            return GF_REACTION_PARSING_ERROR;
+        } catch (const gridfire::exceptions::ReactionError& e) {
+            ctx->last_error = e.what();
+            return GF_REACTION_ERROR;
+        } catch (const gridfire::exceptions::SingularJacobianError& e) {
+            ctx->last_error = e.what();
+            return GF_SINGULAR_JACOBIAN_ERROR;
+        } catch (const gridfire::exceptions::IllConditionedJacobianError& e) {
+            ctx->last_error = e.what();
+            return GF_ILL_CONDITIONED_JACOBIAN_ERROR;
+        } catch (const gridfire::exceptions::CVODESolverFailureError& e) {
+            ctx->last_error = e.what();
+            return GF_CVODE_SOLVER_FAILURE_ERROR;
+        } catch (const gridfire::exceptions::KINSolSolverFailureError& e) {
+            ctx->last_error = e.what();
+            return GF_KINSOL_SOLVER_FAILURE_ERROR;
+        } catch (const gridfire::exceptions::SUNDIALSError& e) {
+            ctx->last_error = e.what();
+            return GF_SUNDIALS_ERROR;
+        } catch (const gridfire::exceptions::SolverError& e) {
+            ctx->last_error = e.what();
+            return GF_SOLVER_ERROR;
+        } catch (const gridfire::exceptions::HashingError& e) {
+            ctx->last_error = e.what();
+            return GF_HASHING_ERROR;
+        } catch (const gridfire::exceptions::UtilityError& e) {
+            ctx->last_error = e.what();
+            return GF_UTILITY_ERROR;
+        } catch (const gridfire::exceptions::DebugException& e) {
+            ctx->last_error = e.what();
+            return GF_DEBUG_ERROR;
+        } catch (const gridfire::exceptions::GridFireError& e) {
+            ctx->last_error = e.what();
+            return GF_GRIDFIRE_ERROR;
+        } catch (std::exception& e) {
+            ctx->last_error = e.what();
+            return GF_NON_GRIDFIRE_ERROR;
+        } catch (...) {
+            ctx->last_error = "Unknown error occurred during evolution.";
+            return GF_UNKNOWN_ERROR;
+        }
+    }
+
+    char* gf_get_last_error_message(void* ptr) {
+        const auto* ctx = static_cast<GridFireContext*>(ptr);
+        return const_cast<char*>(ctx->last_error.c_str());
+    }
+
+    char* gf_error_code_to_string(const int error_code) {
+        switch (error_code) {
+            case GF_SUCCESS:
+                return const_cast<char*>("GF_SUCCESS");
+            case GF_UNKNOWN_ERROR:
+                return const_cast<char*>("GF_UNKNOWN_ERROR");
+            case GF_NON_GRIDFIRE_ERROR:
+                return const_cast<char*>("GF_NON_GRIDFIRE_ERROR");
+            case GF_INVALID_QSE_SOLUTION_ERROR:
+                return const_cast<char*>("GF_INVALID_QSE_SOLUTION_ERROR");
+            case GF_FAILED_TO_PARTITION_ENGINE_ERROR:
+                return const_cast<char*>("GF_FAILED_TO_PARTITION_ENGINE_ERROR");
+            case GF_NETWORK_RESIZED_ERROR:
+                return const_cast<char*>("GF_NETWORK_RESIZED_ERROR");
+            case GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR:
+                return const_cast<char*>("GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR");
+            case GF_BAD_COLLECTION_ERROR:
+                return const_cast<char*>("GF_BAD_COLLECTION_ERROR");
+            case GF_BAD_RHS_ENGINE_ERROR:
+                return const_cast<char*>("GF_BAD_RHS_ENGINE_ERROR");
+            case GF_STALE_JACOBIAN_ERROR:
+                return const_cast<char*>("GF_STALE_JACOBIAN_ERROR");
+            case GF_UNINITIALIZED_JACOBIAN_ERROR:
+                return const_cast<char*>("GF_UNINITIALIZED_JACOBIAN_ERROR");
+            case GF_UNKNOWN_JACOBIAN_ERROR:
+                return const_cast<char*>("GF_UNKNOWN_JACOBIAN_ERROR");
+            case GF_JACOBIAN_ERROR:
+                return const_cast<char*>("GF_JACOBIAN_ERROR");
+            case GF_ENGINE_ERROR:
+                return const_cast<char*>("GF_ENGINE_ERROR");
+            case GF_MISSING_BASE_REACTION_ERROR:
+                return const_cast<char*>("GF_MISSING_BASE_REACTION_ERROR");
+            case GF_MISSING_SEED_SPECIES_ERROR:
+                return const_cast<char*>("GF_MISSING_SEED_SPECIES_ERROR");
+            case GF_MISSING_KEY_REACTION_ERROR:
+                return const_cast<char*>("GF_MISSING_KEY_REACTION_ERROR");
+            case GF_POLICY_ERROR:
+                return const_cast<char*>("GF_POLICY_ERROR");
+            case GF_REACTION_PARSING_ERROR:
+                return const_cast<char*>("GF_REACTION_PARSING_ERROR");
+            case GF_REACTION_ERROR:
+                return const_cast<char*>("GF_REACTION_ERROR");
+            case GF_SINGULAR_JACOBIAN_ERROR:
+                return const_cast<char*>("GF_SINGULAR_JACOBIAN_ERROR");
+            case GF_ILL_CONDITIONED_JACOBIAN_ERROR:
+                return const_cast<char*>("GF_ILL_CONDITIONED_JACOBIAN_ERROR");
+            case GF_CVODE_SOLVER_FAILURE_ERROR:
+                return const_cast<char*>("GF_CVODE_SOLVER_FAILURE_ERROR");
+            case GF_KINSOL_SOLVER_FAILURE_ERROR:
+                return const_cast<char*>("GF_KINSOL_SOLVER_FAILURE_ERROR");
+            case GF_SUNDIALS_ERROR:
+                return const_cast<char*>("GF_SUNDIALS_ERROR");
+            case GF_SOLVER_ERROR:
+                return const_cast<char*>("GF_SOLVER_ERROR");
+            case GF_HASHING_ERROR:
+                return const_cast<char*>("GF_HASHING_ERROR");
+            case GF_UTILITY_ERROR:
+                return const_cast<char*>("GF_UTILITY_ERROR");
+            case GF_DEBUG_ERROR:
+                return const_cast<char*>("GF_DEBUG_ERROR");
+            case GF_GRIDFIRE_ERROR:
+                return const_cast<char*>("GF_GRIDFIRE_ERROR");
+            case FDSSE_NON_4DSTAR_ERROR:
+                return const_cast<char*>("FDSSE_NON_4DSTAR_ERROR");
+            case FDSSE_UNKNOWN_ERROR:
+                return const_cast<char*>("FDSSE_UNKNOWN_ERROR");
+            case FDSSE_SUCCESS:
+                return const_cast<char*>("FDSSE_SUCCESS");
+            case FDSSE_UNKNOWN_SYMBOL_ERROR:
+                return const_cast<char*>("FDSSE_UNKNOWN_SYMBOL_ERROR");
+            case FDSSE_SPECIES_ERROR:
+                return const_cast<char*>("FDSSE_SPECIES_ERROR");
+            case FDSSE_INVALID_COMPOSITION_ERROR:
+                return const_cast<char*>("FDSSE_INVALID_COMPOSITION_ERROR");
+            case FDSSE_COMPOSITION_ERROR:
+                return const_cast<char*>("FDSSE_COMPOSITION_ERROR");
+            default:
+                return const_cast<char*>("GF_UNRECOGNIZED_ERROR_CODE");
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/extern/meson.build b/src/extern/meson.build
new file mode 100644
index 00000000..ba77e65a
--- /dev/null
+++ b/src/extern/meson.build
@@ -0,0 +1,24 @@
+gridfire_extern_sources = files(
+    'lib/gridfire_context.cpp',
+    'lib/gridfire_extern.cpp',
+)
+
+gridfire_extern_dependencies = [
+    gridfire_dep
+]
+
+libgridfire_extern = library('gridfire_extern',
+    gridfire_extern_sources,
+    include_directories: include_directories('include'),
+    dependencies: gridfire_extern_dependencies,
+    install : true
+)
+
+gridfire_extern_dep = declare_dependency(
+    include_directories: include_directories('include'),
+    link_with: libgridfire_extern,
+    sources: gridfire_extern_sources,
+    dependencies: gridfire_extern_dependencies,
+)
+
+install_subdir('include/gridfire', install_dir: get_option('includedir'))
diff --git a/src/meson.build b/src/meson.build
index 89750924..bf6b7014 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -60,6 +60,8 @@ gridfire_dep = declare_dependency(
 
 install_subdir('include/gridfire', install_dir: get_option('includedir'))
 
+subdir('extern')
+
 if get_option('build-python')
     message('Configuring Python bindings...')
     subdir('python')
diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap
index 7048ec12..b97019bc 100644
--- a/subprojects/fourdst.wrap
+++ b/subprojects/fourdst.wrap
@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/fourdst
-revision = v0.9.6
+revision = v0.9.7
 depth = 1
diff --git a/tests/extern/C/gridfire_evolve.c b/tests/extern/C/gridfire_evolve.c
new file mode 100644
index 00000000..c171e7f2
--- /dev/null
+++ b/tests/extern/C/gridfire_evolve.c
@@ -0,0 +1,71 @@
+#include "gridfire/extern/gridfire_extern.h"
+#include <stdio.h>
+#define NUM_SPECIES 8
+
+// Define a macro to check return codes
+#define CHECK_RET_CODE(ret, ctx, msg) \
+    if (ret != 0 && ret != 1) { \
+        printf("Error %s: %s\n", msg, gf_get_last_error_message(ctx)); \
+        gf_free(ctx); \
+        return 1; \
+    }
+
+int main() {
+    void* ctx = gf_init();
+
+    const char* species_names[NUM_SPECIES];
+    species_names[0] = "H-1";
+    species_names[1] = "He-3";
+    species_names[2] = "He-4";
+    species_names[3] = "C-12";
+    species_names[4] = "N-14";
+    species_names[5] = "O-16";
+    species_names[6] = "Ne-20";
+    species_names[7] = "Mg-24";
+    const double abundances[NUM_SPECIES] = {0.702616602672027, 9.74791583949078e-06, 0.06895512307276903, 0.00025, 7.855418029399437e-05, 0.0006014411598306529, 8.103062886768109e-05, 2.151340851063217e-05};
+
+    int ret = gf_register_species(ctx, NUM_SPECIES, species_names);
+    CHECK_RET_CODE(ret, ctx, "SPECIES");
+
+    ret = gf_construct_engine_from_policy(ctx, "MAIN_SEQUENCE_POLICY", abundances, NUM_SPECIES);
+    CHECK_RET_CODE(ret, ctx, "MAIN_SEQUENCE_POLICY");
+
+    ret = gf_construct_solver_from_engine(ctx, "CVODE");
+    CHECK_RET_CODE(ret, ctx, "CVODE");
+
+    double Y_out[NUM_SPECIES];
+    double energy_out;
+    double dEps_dT;
+    double dEps_dRho;
+    double mass_lost;
+
+    ret = gf_evolve(
+        ctx,
+        abundances,
+        NUM_SPECIES,
+        1.5e7,    // Temperature in K
+        1.5e2,    // Density in g/cm^3
+        3e17,      // Time step in seconds
+        1e-12, // Initial time step in seconds
+        Y_out,
+        &energy_out,
+        &dEps_dT,
+        &dEps_dRho, &mass_lost
+    );
+
+    CHECK_RET_CODE(ret, ctx, "EVOLUTION");
+
+
+    printf("Evolved abundances:\n");
+    for (size_t i = 0; i < NUM_SPECIES; i++) {
+        printf("Species %s: %e\n", species_names[i], Y_out[i]);
+    }
+    printf("Energy output: %e\n", energy_out);
+    printf("dEps/dT: %e\n", dEps_dT);
+    printf("dEps/dRho: %e\n", dEps_dRho);
+    printf("Mass lost: %e\n", mass_lost);
+
+    gf_free(ctx);
+
+    return 0;
+}
\ No newline at end of file
diff --git a/tests/extern/C/meson.build b/tests/extern/C/meson.build
new file mode 100644
index 00000000..80db9f58
--- /dev/null
+++ b/tests/extern/C/meson.build
@@ -0,0 +1,5 @@
+executable('test_c_extern', 'gridfire_evolve.c',
+  install: false,
+  c_args: ['-Wall', '-Wextra'],
+  dependencies: [gridfire_extern_dep]
+)
\ No newline at end of file
diff --git a/tests/extern/meson.build b/tests/extern/meson.build
new file mode 100644
index 00000000..ed6c93a5
--- /dev/null
+++ b/tests/extern/meson.build
@@ -0,0 +1 @@
+subdir('C')
\ No newline at end of file
diff --git a/tests/meson.build b/tests/meson.build
index 6914a261..a2b8027c 100644
--- a/tests/meson.build
+++ b/tests/meson.build
@@ -5,3 +5,4 @@ gtest_nomain_dep = dependency('gtest', main: false, required : true)
 
 # Subdirectories for unit and integration tests
 subdir('graphnet_sandbox')
+subdir('extern')