guidance on code structure improvement

I would like to know if this code can be improved/optimized. is it well-written? Any help will be appreciated. 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

[dpb]

Comments:

1. No comments -- you may remember what you did/are doing tomorrow, but six months down the road...or the person who may have to pick up after you?
2. I'd remove the many extra blank lines -- they just take up room for seeing the code flow succinctly. (Some seem to like a lot of white space, I don't).
3. I'd use a Switch block here for a selection of method instead of the if...else...end
4. Would be a more generic function and wouldn't require changing code to pass the inputs to the top level function. You could assign default values if user doesn't pass some/any.
5. What's the point of the do-nothing expressions 'y+0' in function solve_radau_adaptive? It also has a bunch of unused parameters in the argument list -- maybe it's not complete?