Intuitive GSL stiff EDO

Question

I need to solve an ODE that has stiff behaviour. Said ODE would be given by

int odefunc (double x, const double y[], double f[], void *params)
{
    double k=1e12; //At k=1e7 it works (doesn't show stiffness)
    double y_eq=-1.931+1.5*log(x)-x;
    f[0] = k*(exp(2*y_eq-y[0])-exp(y[0]));   
    return GSL_SUCCESS;
}

I'm not so used to C (I am to C++, though) and the documentation in the GSL page for Differential Equations doesn't have that many comments for me to understand it. I've used this post to understand it (also because I only have 1 ode, not a system of them), and to some extent it has been useful, but it's still a bit confusing to me. For my tests I've used my function instead of the one given, and I've gone with it although I don't understand the code. But now it was key for me to understand, because...

I would need to modify the examples they give me, I need to use an algorithm that solves stiff equations, preferably gsl_odeiv2_step_msbdf with an adaptative step-size or something equivalent to that (the BDF method seems to be used profusely in the academia). That needs the jacobian, and that part is really intricate if you're not used to GSL or C.

Answer 1

To manually implement an algorithmic differentiation it helps to include explicit intermediary steps

int odefunc (double x, const double y[], double f[], void *params)
{
    double k=1e12; //At k=1e7 it works (doesn't show stiffness)
    double y_eq=-1.931+1.5*log(x)-x;
    double v1 = exp(2*y_eq-y[0]);
    double v2 = exp(y[0]);
    double z = k*(v1-v2);   
    f[0] = z;
    return GSL_SUCCESS;
}

Then you can enhance each step with its derivative propagation

int jac (double x, const double y[], double *dfdy,
     double dfdx[], void *params)
{
    double k=1e12; //At k=1e7 it works (doesn't show stiffness)
    // Dx_dx=Dy_dy=1, Dx_dy=Dy_dx=0 are used without naming them
    double y_eq = -1.931+1.5*log(x)-x, 
          Dy_eq_dx=1.5/x-1, 
          Dy_eq_dy=0;
    double v1 = exp(2*y_eq-y[0]),
          Dv1_dx=v1*(2*Dy_eq_dx-0),
          Dv1_dy=v1*(2*Dy_eq_dy-1);
    double v2 = exp(y[0]), 
          Dv2_dx=v2*(0), 
          Dv2_dy=v2*(1);
    double z = k*(v1-v2), 
          Dz_dx=k*(Dv1_dx-Dv2_dx), 
          Dz_dy=k*(Dv1_dy-Dv2_dy);   
    dfdx[0] = Dz_dx;
    dfdy[0] = Dz_dy;
    return GSL_SUCCESS;
}

For larger code use a code-rewriting tool that does these steps automatically. auto-diff.org should have a list of appropriate ones.