Implementation of TR-BDF2 in matlab
<p>I want to implement TR-BDF2 adaptive-timestep version in matlab.
The C program is in <a href="https://math.la.asu.edu/~gardner/TRBDF2.c" rel="nofollow noreferrer">https://math.la.asu.edu/~gardner/TRBDF2.c</a> and there are a lot of online C++ version. But for matlab, I cannot find any except the constant-timestep version and thus have difficulty translating the C program to matlab version. 
Also, I hope to implement it in multivariable version in vector-matrix form. What should I pay attention to?</p>

<p>The essence of <a href="https://math.la.asu.edu/~gardner/TRBDF.pdf" rel="nofollow noreferrer">adaptive timestepping</a> is that I have to select each timestep inside a window:</p>

<p><a href="https://i.stack.imgur.com/Hf8fJ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Hf8fJ.png" alt="enter image description here"></a></p>

<p>In C program above it seems that this does it->The TR-BDF2 timestep selection is: </p>

<pre><code>void compute_dt(GRID *grid, TIMESTEP *timestep, double F[], double F_mid[],
    double F_new[])
{
    int modes = grid-&gt;modes, i;
    double e_l; // = ||e_l||
    double dt = timestep-&gt;dt, dt_new;
    double dt_max = timestep-&gt;dt_max;
    double r, k = (-3.*sq(GAMMA)+4.*GAMMA-2.)/(12.*(2.-GAMMA));
    double EPSILON_REL = 1.e-6, EPSILON_ABS = 1.e-12;

    e_l = 0.;
    for (i = 0; i &lt; modes; i++)
        e_l += fabs(2.*k*dt*(F[i]/GAMMA - F_mid[i]/(GAMMA*(1.-GAMMA)) +
            F_new[i]/(1.-GAMMA)));
    e_l /= modes;
    r = e_l/(EPSILON_REL*timestep-&gt;norm_u + EPSILON_ABS);
    if (DEBUG_dt) printf("TRBDF2 r = %g\n", r);

    if (r &lt;= 2.) {
        timestep-&gt;REDO_STEP = NO;
        dt_new = dt/pow(r, 1./3.);
        dt = min(dt_new, 5.*dt);
    }
    else {
        timestep-&gt;REDO_STEP = YES;
        dt /= 2.;
        printf("WARNING: TRBDF2 r = %g &gt; 2\n", r);
        printf("dt /= 2.; redoing timestep %d\n", timestep-&gt;step);
    }
    dt = min(dt, dt_max);

    timestep-&gt;dt = dt;
    timestep-&gt;FACTOR = dt/timestep-&gt;dt_euler;
}

void compute_first_dt(GRID *grid, TIMESTEP *timestep, double u[])
{
    int modes = grid-&gt;modes;
    double *D;

    D = (double *) alloc_vector(modes, DOUBLE);
    compute_D(grid, u, D);
    compute_dt_euler(grid, timestep, D);
    free((char *) D);
    timestep-&gt;dt = timestep-&gt;FACTOR*timestep-&gt;dt_euler;
    timestep-&gt;dt_min = timestep-&gt;MIN_FACTOR*timestep-&gt;dt_euler;
    timestep-&gt;dt_max = timestep-&gt;MAX_FACTOR*timestep-&gt;dt_euler;
}
</code></pre>

<p>Also, this part also do the comparison:</p>

<pre><code>void select_dt(TIMESTEP *timestep)
{
    double t_max = timestep-&gt;t_max;

    timestep-&gt;LAST_STEP = NO;
    if (timestep-&gt;t+timestep-&gt;dt &gt; t_max) {
        timestep-&gt;dt = t_max - timestep-&gt;t;
        timestep-&gt;t = t_max;
        timestep-&gt;FACTOR = timestep-&gt;dt/timestep-&gt;dt_euler;
        timestep-&gt;LAST_STEP = YES;
    }
    else timestep-&gt;t += timestep-&gt;dt;
}
</code></pre>

<p>But how can I do this in matlab? And does this scheme change in multivariable version of ode?</p>