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 ?

Answer 1

[Umar]

Hi @Cesar,

Thanks for sharing your code — it is clear you have a solid grasp of the method and the structure is easy to follow. That said, I went through it carefully and found a few things worth fixing before you take this further.

The most important issue is that your b4 and b5 labels are swapped. In the published Dormand-Prince tableau, b5 should hold the coefficients [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0] and b4 the others — but in your code it is the other way around. The numerical values are all correct, so your error estimate still works, but the line u = y5 ends up advancing the integration with the lower-accuracy 4th-order solution instead of the 5th-order one. Dormand and Prince specifically designed the method so the 5th-order solution propagates forward — that is the whole point of DOPRI over the earlier Fehlberg method — so this is worth fixing first.

The second issue is your error criterion. You are using norm(y5 - y4, inf) <= tol which is a flat absolute threshold. MATLAB's ode45 uses a mixed relative and absolute formula: sc = AbsTol + RelTol * max(abs(y4), abs(y5)), then normalises the error against sc. In the Robertson problem, y2 drops to around 1e-5, so with your tol = 1e-2 the solver is essentially blind to inaccuracies in that component. This is exactly what Torsten is getting at when he asks you to compare step counts against ode45 — the two solvers would diverge immediately because they are not measuring the same thing when deciding whether to accept a step.

On that note, Robertson is actually the wrong problem to validate this solver on. It is a canonically stiff system — MathWorks uses it as a benchmark for ode15s, not ode45. On your short tspan the step size is being squeezed by stability constraints, not accuracy, so you are not really testing your error control at all. A much better validation test, which is what Torsten is suggesting, is to run your solver against ode45 on something simple and non-stiff like y' = -y or a harmonic oscillator, with equal tolerances, and check that your step counts are in the same ballpark. If they are not, Issues 1 and 2 above are the reason.

A few more things worth tidying up:

Your step size growth is using fixed absolute bounds (h_min and h_max) rather than the relative growth cap that ode45 uses. MATLAB's formula ensures the step never grows by more than a factor of roughly 5 or shrinks below 10% in a single update — your current approach does not have that guardrail.

The FSAL property is not implemented. Dormand-Prince was designed so that the last stage of one step (k7) is reused as the first stage of the next, saving one function evaluation per accepted step — about 14% less work overall. Right now you are recomputing k(:,1) = f(t, u) fresh every iteration.

The Robertson system satisfies y1 + y2 + y3 = 1 at all times. It is worth adding a quick check after each accepted step to make sure sum(u) stays close to 1 — it is a free sanity check and with your current loose error control it could be flagging something.

Finally, your error subplot label says the quantity is divided by h but the code never actually divides by h — either fix the label or fix the plot. Also, Figures 2 and 3 are plotting the same data with only minor cosmetic differences, so those can be merged into one.

To summarise the priority order: fix the coefficient swap first, then the error normalisation, then run Torsten's benchmark on a non-stiff problem. Everything else is secondary to those three.

The foundation is genuinely solid — the Butcher tableau values are correct and the adaptive loop logic is coherent. These are targeted fixes, not a rewrite.

Good luck with it.

[Cesar Cardenas]

Thank you for your detailed response. I have this code now, can it be improved, not sure if it'es well-written? Thanks

function main_solver_optimized()

    params.alpha = 0.04; params.beta = 1e4; params.gamma = 3e7;
    y0 = [1; 0; 0];
    tspan = [0, 1e6];
    h0 = 1e-4;
    tol = 1e-2;

 
    fprintf('Running Part 1: Radau IIA (2-stage vs 3-stage)...\n');
    [T1, Y1, s1, r1] = solve_robertson(y0, tspan, h0, tol, params, 'Radau');

    
    fprintf('Running Part 2: SDIRK 3(4)...\n');
    [T2, Y2, s2, r2] = solve_robertson(y0, tspan, h0, tol, params, 'SDIRK');


    fprintf('\n================ RESULTS COMPARISON ================\n');
    fprintf('Method       | Successful Steps | Rejected Steps\n');
    fprintf('----------------------------------------------------\n');
    fprintf('Radau IIA    | %16d | %14d\n', s1, r1);
    fprintf('SDIRK 3(4)   | %16d | %14d\n', s2, r2);
    fprintf('====================================================\n');

  
    figure('Name', 'Robertson Problem Solutions');
    subplot(2,1,1);
    semilogx(T1, Y1(:,1), T1, Y1(:,2)*1e4, T1, Y1(:,3), 'LineWidth', 1.5);
    title('Part 1: Radau IIA Method'); xlabel('Time (s)'); ylabel('Concentration');
    legend('y_1', 'y_2 \times 10^4', 'y_3', 'Location', 'Best'); grid on;

    subplot(2,1,2);
    semilogx(T2, Y2(:,1), T2, Y2(:,2)*1e4, T2, Y2(:,3), 'LineWidth', 1.5);
    title('Part 2: SDIRK 3(4) Method'); xlabel('Time (s)'); ylabel('Concentration');
    legend('y_1', 'y_2 \times 10^4', 'y_3', 'Location', 'Best'); grid on;
end


function [T, Y, success, rejected] = solve_robertson(y0, tspan, h, tol, p, method)
   
    Nmax = 1e6;
    T = zeros(Nmax,1); Y = zeros(Nmax,3);
    t = tspan(1); y = y0;
    T(1) = t; Y(1,:) = y';
    stepCount = 1;
    success = 0; rejected = 0;


    if strcmp(method,'Radau')
   
        A2 = [5/12, -1/12; 3/4, 1/4]; c2 = [1/3, 1]';
       
        sq6 = sqrt(6);
        A3 = [(88-7*sq6)/360, (296-169*sq6)/1800, (-2+3*sq6)/225; ...
              (296+169*sq6)/1800, (88+7*sq6)/360, (-2-3*sq6)/225; ...
              (16-sq6)/36, (16+sq6)/36, 1/9];
        c3 = [(4-sq6)/10, (4+sq6)/10, 1]';
        solver_func = @(y,h) solve_radau_adaptive(y,h,A2,c2,A3,c3,p);
        order_p = 3;
    else
    
        gam = 1/4;
        A = [gam, 0, 0, 0, 0; 
             1/2-gam, gam, 0, 0, 0; 
             2*gam, 1-4*gam, gam, 0, 0; 
             1/10, 1/10, 1/10+gam, gam, 0; 
             1/4, 1/4, 0, 1/4, gam];
        c = [1/4, 3/4, 11/20, 1/2, 1]';
        b = [1/4, 1/4, 0, 1/4, 1/4];      
        bhat = [-3/16, 5/8, 3/16, 1/16, 5/16]; 
        solver_func = @(y,h) solve_sdirk_sequential(y,h,A,b,bhat,c,p);
        order_p = 3;
    end

    
    while t < tspan(2)
        if t + h > tspan(2), h = tspan(2) - t; end
        [y_n, y_hat, converged] = solver_func(y,h);

        if converged
            err = max(abs(y_n - y_hat));
            if err <= tol
                t = t + h; y = y_n;
                stepCount = stepCount + 1;
                T(stepCount) = t;
                Y(stepCount,:) = y';
                success = success + 1;

               
                h = h * min(5, max(0.1, 0.9*(tol/err)^(1/(order_p+1))));
            else
                h = h/4; rejected = rejected + 1;
            end
        else
            h = h/4; rejected = rejected + 1;
        end
    end

  
    T = T(1:stepCount); Y = Y(1:stepCount,:);
end

function [ynext, yhat, conv] = solve_radau_adaptive(y, h, A2, c2, A3, c3, p)

    ynext = y + 0;
    yhat  = y + 0;
    conv = true;
   
end

function [ynext, yhat, conv] = solve_sdirk_sequential(y, h, A, b, bhat, c, p)
    s = length(c);
    K = zeros(3,s);
    conv = true;
    for i = 1:s
        ki = zeros(3,1); 
        rhs = y + h*(K(:,1:i-1)*A(i,1:i-1)');
        for iter = 1:5
            val = rhs + h*A(i,i)*ki;
            f_val = [-p.alpha*val(1) + p.beta*val(2)*val(3); ...
                      p.alpha*val(1) - p.beta*val(2)*val(3) - p.gamma*val(2)^2; ...
                      p.gamma*val(2)^2];
            jac = [-p.alpha, p.beta*val(3), p.beta*val(2); ...
                    p.alpha, -p.beta*val(3)-2*p.gamma*val(2), -p.beta*val(2); ...
                    0, 2*p.gamma*val(2), 0];
            
            F = ki - f_val;
            DF = eye(3) - h*A(i,i)*jac;
            step = DF\F;
            ki = ki - step;
            if max(abs(step)) < 1e-8, break; end
            if iter == 5, conv = false; end
        end
        K(:,i) = ki;
    end
    ynext = y + h*(K*b');
    yhat  = y + h*(K*bhat');
end

[Umar]

Hi @Cesar,

Good progress overall — moving to Radau IIA and SDIRK was the right call. Robertson is a stiff problem and implicit methods are exactly what it needs. The Butcher tableau values are correct and the SDIRK Newton solver has solid structure.

That said, there are two things stopping this from producing real results right now.

The Radau solver doesn't actually do anything. It just returns the initial condition on every step. So the entire Part 1 comparison is meaningless — you're printing step counts for a solver that never integrated a single point.

The error criterion is still the same flat threshold from your last version. The Robertson y2 component drops to around 1e-5, and you're tolerating errors up to 1e-2. The solver is basically blind to what's happening in that component.

Fix those two things first and everything else becomes secondary. The Newton loop also needs a proper early break rather than running all 5 iterations before flagging failure, and it's worth double-checking the bhat coefficients against Hairer & Wanner since a negative weight on the first stage can cause the error estimator to overreact on stiff problems.

The foundation is there — you just can't trust any of the numbers it's printing yet.