code improvement/optimization

I would like to know if this code can be improved or it well written? any comments will be greatly appreciated. Thanks

u0 = [1; 0; 0];
tspan = [0, 1e-3];   % small final time for explicit solver
tol = 1e-2;
h0 = 1e-6;

% Execute solver
[t, u, err, h_steps, stats] = dormand_prince_solver(@robertson_ode, tspan, u0, tol, h0);

% Part (a): Plot Solutions
figure;
subplot(2,1,1);
plot(t, u(:,1), 'b', t, u(:,3), 'g', 'LineWidth', 1.5);
legend('y1 (Reactant)', 'y3 (Product)');
xlabel('Time t'); ylabel('Concentration');
%title('Robertson Problem: Concentrations y1 and y3'); 
grid on;

subplot(2,1,2);
semilogx(t, u(:,2), 'r', 'LineWidth', 1.5);
xlabel('Time t (Log Scale)'); ylabel('Concentration y2');
%title('Robertson Problem: Intermediate Species y2'); 
grid on;

% Part (b): Error and Stepsize Plots
figure;
subplot(2,1,1);
semilogy(t(2:end), err);
xlabel('Time t'); ylabel('||y_n^5 - y_n^4||_{\infty} / h');
%title('Part (b): Estimated Error vs Time'); 
grid on;

subplot(2,1,2);
semilogy(t(2:end), h_steps);
xlabel('Time t'); ylabel('Stepsize h');
%title('Part (b): Stepsize h vs Time'); 
grid on;

figure;

subplot(2,1,1);
semilogy(t(2:end), err, 'LineWidth', 1.5);
xlabel('Time t');
ylabel('Error');
title('Estimated Local Error vs Time');
grid on;

subplot(2,1,2);
semilogy(t(2:end), h_steps, 'LineWidth', 1.5);
xlabel('Time t');
ylabel('Stepsize h');
title('Stepsize vs Time');
grid on;

function dydt = robertson_ode(~, y)
    % Coefficients from Screenshot 2026-02-28 at 3.38.10 PM.png
    alpha = 0.04;
    beta = 1e4;
    gamma = 3e7;
    
    dydt = [ -alpha*y(1) + beta*y(2)*y(3);
              alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)^2;
              gamma*y(2)^2 ];
end

function [t_all, u_all, error_history, h_history, stats] = ...
    dormand_prince_solver(f, tspan, u0, tol, h0)

    % Initialization
    t = tspan(1);
    u = u0(:);
    h = h0;

    t_all = t;
    u_all = u.';
    error_history = [];
    h_history = [];

    n_accepted = 0;
    n_rejected = 0;

    % ---- Dormand–Prince 4(5) Coefficients ----
    a = [1/5,         0,          0,          0,         0,          0;
         3/40,        9/40,       0,          0,         0,          0;
         44/45,      -56/15,      32/9,       0,         0,          0;
         19372/6561, -25360/2187, 64448/6561, -212/729,  0,          0;
         9017/3168,  -355/33,     46732/5247, 49/176,    -5103/18656, 0;
         35/384,      0,          500/1113,   125/192,   -2187/6784,  11/84];

    b4 = [5179/57600, 0, 7571/16695, 393/640, ...
         -92097/339200, 187/2100, 1/40];

    b5 = [35/384, 0, 500/1113, 125/192, ...
         -2187/6784, 11/84, 0];

    c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];

    safety = 0.9;
    h_min = 1e-12;
    h_max = 1.0;

    while t < tspan(2)

        if h < h_min
            error('Step size underflow: problem likely stiff.');
        end

        if t + h > tspan(2)
            h = tspan(2) - t;
        end

        % ---- Stage calculations ----
        k = zeros(length(u), 7);
        k(:,1) = f(t, u);

        for i = 2:7
            k(:,i) = f(t + c(i)*h, ...
                u + h * (k(:,1:i-1) * a(i-1,1:i-1)'));
        end

        % ---- 4th and 5th order solutions ----
        y4 = u + h * (k * b4');
        y5 = u + h * (k * b5');

 
        error_val = norm(y5 - y4, inf);

   
        if error_val <= tol

     
            t = t + h;
            u = y5;

            t_all = [t_all; t];
            u_all = [u_all; u.'];
            error_history = [error_history; error_val];
            h_history = [h_history; h];

            n_accepted = n_accepted + 1;

        else
            n_rejected = n_rejected + 1;
        end

      
        if error_val == 0
            h_new = 2*h;
        else
            h_new = h * safety * (tol / error_val)^(1/5);
        end

        h = min(max(h_new, h_min), h_max);
    end

    stats.accepted = n_accepted;
    stats.rejected = n_rejected;
end

[Torsten]

Did you test the code for a non-stiff system of differential equations against ode45 for a reasonable timespan ? Are the number of steps and the results similar for equally chosen tolerances ?