odexprb34 solver

As I did with the odexprk23 solver from the Runge-Kutta family, I implemented the official exponential fourth order method with adaptive time stepping from Rosenbrock family described in my post.  Also in this case I used the function expmv and I wrote a code for the stepper (exp_rosenbrock_34), 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 + \frac{1}{2} h_n \varphi_1( \frac{1}{2} h_n J_n) F(t_n,u_n) + \frac{1}{2} h_n \varphi_2(\frac{1}{2} h_n J_n) \frac{1}{2} h_n v_n ,
\\
U_{n3} = u_n +  h_n \varphi_1( h_n J_n) [F(t_n,u_n) + D_{n2}] + h_n \varphi_2 ( h_n J_n) h_n v_n ,
$$
[Atilde, eta] = (1/2*dt, Jn, [Fn, 1/2*dt*vn]);
X = myexpmv (1/2*dt, Atilde, [zeros(size (Atilde)-1, 1); 1/eta])(1:d, :);
U(:, 2) = x(:) + X;

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

The third call to expmv is for 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 J_n) F(t_n,u_n) + h_n \varphi_2(h_n J_n) h_n v_n \dots
\\
\dots + h_n \varphi_3(h_n J_n) [ 16 D_{n2} - 2 D_{n3}] + h_n \varphi_4(h_n J_n) [ -48 D_{n2} + 12 D_{n3}] ,
\\
\hat{u}_{n+1} =  u_n + h_n \varphi_1(h_n J_n) F(t_n,u_n) + h_n \varphi_2(h_n J_n) h_n v_n + h_n \varphi_3(h_n J_n) [ 16 D_{n2} - 2 D_{n3}] ,
$$
v = [Fn, dt*vn, 16*D(:, 2)-2*D(:, 3), -48*D(:, 2)+12*D(:, 3)];
w = [Fn, dt*vn, 16*D(:, 2)-2*D(:, 3)];
[Atilde, eta] = augmat (dt, Jn, 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 (odexprb34), according to odepkg standards and  I tested the order of the method with the usual example described in my post.