1. Wrap the necessary headers to sundials.c:

#define _Float16 float // don't care why we don't find _Float16 on Apple Silicon

#include "cvode/cvode.h"
#include "nvector/nvector_serial.h"
#include "sunmatrix/sunmatrix_dense.h"
#include "sunlinsol/sunlinsol_dense.h"

2. Port the cvRoberts_dns.c example to D:

import sundials;
import std.stdio;

extern (C) static int f(sunrealtype t,
                        N_Vector y,
                        N_Vector ydot,
                        void* user_data)
{
  sunrealtype y1, y2, y3, yd1, yd3;

  y1 = NV_DATA_S(y)[0];
  y2 = NV_DATA_S(y)[1];
  y3 = NV_DATA_S(y)[2];

  yd1 = NV_DATA_S(ydot)[0] = -0.04 * y1 + 1.0e4 * y2 * y3;
  yd3 = NV_DATA_S(ydot)[2] = 3.0e7 * y2 * y2;
  NV_DATA_S(ydot)[1]       = -yd1 - yd3;

  return (0);
}

extern (C) static int jac(sunrealtype t,
                          N_Vector y,
                          N_Vector fy,
                          SUNMatrix J,
                          void* user_data,
                          N_Vector tmp1,
                          N_Vector tmp2,
                          N_Vector tmp3)
{
    double y2 = NV_DATA_S(y)[1];
    double y3 = NV_DATA_S(y)[2];

    auto J_content = cast(SUNMatrixContent_Dense) J.content;

    J_content.cols[0][0] = -0.04;
    J_content.cols[1][0] = 1.0e4 * y3;
    J_content.cols[2][0] = 1.0e4 * y2;

    J_content.cols[0][1] = 0.04;
    J_content.cols[1][1] = -1.0e4 * y3 - 6.0e7 * y2;
    J_content.cols[2][1] = -1.0e4 * y2;

    J_content.cols[0][2] = 0.0;
    J_content.cols[1][2] = 6.0e7 * y2;
    J_content.cols[2][2] = 0.0;

    return (0);
}

extern(C) static int g(sunrealtype t,
                       N_Vector y,
                       double* gout,
                       void* user_data)
{
    double y1 = NV_DATA_S(y)[0];
    double y3 = NV_DATA_S(y)[2];

    gout[0] = y1 - 0.0001;
    gout[1] = y3 - 0.01;

    return 0;
}

void main()
{
    SUNContext sunctx;
    int retval = SUNContext_Create(SUN_COMM_NULL, &sunctx);
    writeln("SUNContext created with return value ", retval);

    N_Vector y = N_VNew_Serial(3, sunctx);
    NV_DATA_S(y)[0] = 1.0;
    NV_DATA_S(y)[1] = 0.0;
    NV_DATA_S(y)[2] = 0.0;
    writeln("State vector created");

    void* cvode_mem = CVodeCreate(CV_BDF, sunctx);
    retval = CVodeInit(cvode_mem, &f, 0.0, y);
    writeln("cvode_mem initialized, return value is ", retval);


    SUNMatrix A = SUNDenseMatrix(3, 3, sunctx);

    SUNLinearSolver LS = SUNLinSol_Dense(y, A, sunctx);
    retval = CVodeSetLinearSolver(cvode_mem, LS, A);
    writeln("Linear solver set, return value is ", retval);

    retval = CVodeSetJacFn(cvode_mem, &jac);
    writeln("Jacobian function set, return value is ", retval);

    N_Vector abstol = N_VNew_Serial(3, sunctx);
    NV_DATA_S(abstol)[0] = 1.0e-8;
    NV_DATA_S(abstol)[1] = 1.0e-14;
    NV_DATA_S(abstol)[2] = 1.0e-6;
    writeln("Absolute tolerance vector created");
    retval = CVodeSVtolerances(cvode_mem, 1.0e-4, abstol);

    retval = CVodeRootInit(cvode_mem, 2, &g);
    writeln("Root finding initialized, return value is ", retval);

    double t = 0.0;
    double tout = 0.4;

    int[2] roots_found;
    int n_out = 12;
    int i_out = 0;

    double t_mult = 10.0;

    while(true)
    {
        retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
        writeln(t, "\t", NV_DATA_S(y)[0], "\t", NV_DATA_S(y)[1], "\t", NV_DATA_S(y)[2]);

        if (retval == CV_ROOT_RETURN)
        {
            retval = CVodeGetRootInfo(cvode_mem, roots_found.ptr);
            writeln("Root found at t = ", t, " with g1 = ", roots_found[0], " and g2 = ", roots_found[1]);
        }

        if (retval == CV_SUCCESS)
        {
            i_out++;
            tout *= t_mult;
        }

        if (i_out == n_out) break;
    }
}

(Pls pardon my French, this has been done in cca. 15 minutes)

3. Build:

ldc2 sundials_example.d sundials.c -L-L/usr/local/lib -L-lsundials_core -L-lsundials_nvecserial -L-lsundials_cvode

4. Enjoy!

This. Is. Amazing.

Thank you for all involved in importC development!