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

Answer 1

[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?

Answer 2

[Rik]

The most important thing is documentation. Code without documentation is usefull exactly once: right now, right after you've written it.

For every function you should write a short description of what it does, along with syntax examples and expected output. For each Thing™ in your code you should have a comment. That way you write your functions in a modular way (they can be replaced by a better version without the rest of the code needing any change) and you write your code itself in a modular way (you or someone else can understand what each section/line is doing and can try to fix a bug or implement something faster).

I also think this is too much whitespace. I prefer at most 1 blank line between parts and I don't like to indent the content of my functions. I once followed a course where the instructor insisted on an 8-space indent and a maximum of 80 characters per line. His reason was that it forces you to carefully think which sections you could separate into standalone functions. While I think that is rather extreme (I use 4-space and 100), there is something to be said about the idea.

[dpb]

_" I don't like to indent the content of my functions. . the instructor insisted on an 8-space indent and a maximum of 80 characters per line. ... (I use 4-space and 100),..."

I like to indent the function content. I also typically set line length at 100 but am not excessively reliouious about occasionally allowing a comment to go to about 120. The (one and only <g>) indentation depth is 2 spaces -- enough to be clear is indented yet doesn't waste all that precious real estate at the LHS of the page otherwise. The olden days of 80 column punch cards and text-only 14" monochrome monitors are long past now.

My commenting style other than the preamble as @Rik described is almost all at end of code line(*) and I work diligently to keep the starting point of those aligned over sections of code of roughly equal length lines, not just at the end of the line. A separate section may have a prefix comment line that describes its purpose while a particularly difficult or clever-coding section may have a more detailed synopsis.

(*) And comments must actually add something to the explanation -- a line such as one sees all too frequently x=3; % set x to 3

is useless and is better left off. WHY 3 and who is x? are the important points.

One other stylistic preference since we're on a roll -- use succinct but obvious abbreviations for variable names. LongWindedandexcessivelywordy variables are both error prone for typos and just confuse more than help by making the code exceedingly dense in characters to try to find the operators/logic.

In a similar vein, don't clutter arithmetic expressions with needless parentheses--use them as needed for grouping, of course, but avoid nesting deeper than needed. 2*x.*exp(y)/(1+3z) doesn't need ((2*x).*exp(y))/(1+3z) or somesuch as one frequently sees.

I'm also compulsive in data -- in the allocation of the arrays, I'd align the fields as

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];

as it makes verifying easier. It's somewhat like the mentality of the electrician who dresses hs cable runs in wall cavities.