odexprk23 solver

After the midterm evaluation, I implemented the official exponential method that will be in odepkg, it is the advanced third order method with adaptive time stepping from Runge-Kutta family described in my post. Unlike the general schemes, for this code I used the function expmv so that I could evaluate a single exponential of an augmented matrix instead of computing any φ functions. I wrote a code for the stepper (exp_runge_kutta_23), in which there are three calls to expmv function. The first two for the calculation of the second and third stage
$$
U_{n2} = u_n + c_2 h_n \varphi_1(c_2 h_n A) F(t_n,u_n) ,
\\
U_{n3} = u_n + \frac{2}{3} h_n \varphi_1(\frac{2}{3} h_n A) F(t_n,u_n) + \frac{2}{3} h_n \varphi_2(\frac{2}{3} h_n A) \frac{2}{3 c_2}D_{n2} ,
$$
myexpmv = @(t,A,v) expmv (t, A, v);

[Atilde, eta] = augmat (c2*dt, A, Fn);
X = myexpmv (c2*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);
U(:, 2) = x(:) + X;

[Atilde, eta] = augmat (2/3*dt, A, [Fn, 2/(3*c2)*D(:, 2)]);
X = myexpmv (2/3*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);
U(:,3) = x(:) + X;

Then I changed the file augmat.m to calculate a "double" augmented matrix, and in this way, I have been able to calculate, with a single call to expmv, both the solution to the next step, and the lower order solution for the estimation of the error

$$
u_{n+1} = u_n + h_n \varphi_1(h_n A) F(t_n,u_n) + h_n \varphi_2(h_n A) \frac{3}{2}D_{n3} ,
\\
\hat{u}_{n+1} = u_n + h_n \varphi_1(h_n A) [ F(t_n,u_n) + \frac{1}{2 c_2}D_{n2}],
$$
v = [Fn, 3/2*D(:, 3)];
w = [Fn + 1/(2*c2)*D(:, 2)];
[Atilde, eta] = augmat (dt, A, v, w);
p = size (v, 2);
q = size (w, 2);
ep = [[zeros(p-1, 1); 1/eta; zeros(q, 1)], [zeros(p+q-1, 1); 1/eta]];
X = [eye(d, d), zeros(d, p+q)]*myexpmv(dt, Atilde, [sparse(d, 1), sparse(d, 1); ep]);
X1 = X(:, 1);
X2 = X(:, 2);
x_next = x(:) + X1;
x_est = x(:) + X2;

I wrote also a code for the solver (odexprk23), according to odepkg standards.
I finally tested the order of the method with the usual example described in my post.