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

[Cesar Cardenas]

Thank you for your feedback. Can some of you gvive nay comments on this code? can it also be improved? Thanks

function [t, Y] = mid_term_pc(~, ~, ~, ~)

        B = [-1 0 0; 0 -2 0; 0 0 -3];
        y0 = [1;1;1];
        tspan = [0 5];
        h = 0.1;
  


    t0 = tspan(1);
    tf = tspan(2);
    
    t = t0:h:tf;
    N = length(t);
    
    m = length(y0);
    Y = zeros(m, N);
    Y(:,1) = y0;
    
    A = [5/12, -1/12;
         3/4,   1/4];
    b = [3/4; 1/4];

    I = eye(m);

   
    M = [I - h*A(1,1)*B,   -h*A(1,2)*B;
         -h*A(2,1)*B,      I - h*A(2,2)*B];

    
    [L,U,P] = lu(M);

    for n = 1:N-1
        yn = Y(:,n);

        
        F = [B*yn;
             B*yn];

       
        K = U \ (L \ (P*F));

 
        k1 = K(1:m);
        k2 = K(m+1:end);


        Y(:,n+1) = yn + h*(b(1)*k1 + b(2)*k2);
    end

    % Plot
plot(t, Y)
xlabel('t')
ylabel('y components')
legend('y1','y2','y3')
title('Radau 2-stage solution')
end

[dpb]

I would suggest you apply the previoius comments to the code first...

Answer 1

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

Edit to remove ambiguity: I do use indentation inside functions, but I don't indent the entire function content simply for being in a function. That smells to much like Python to me, and I like to keep my languages separated. I guess it is one of those weird things you develop if you hang around for long enough.

[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 religious about not occasionally allowing a comment to go to about 120 (I turn off the auto code comment wrapping in the editor).. 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 and I can't think of the last time I actually printed a code listing so that no longer is an issue at least for me.

(+) Again my practice here isn't 100% consistent, I may choose to keep the LH real estate rather than indenting the whole function code.

My commenting style other than the preamble as @Rik described is almost all at end of code lines(*) 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.

Answer 2

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

[Cesar Cardenas]

I have added some comments.Thank yiou

function [t, Y] = mid_term_pc(B, y0, tspan, h, doPlot)
% Solution for y' = B*y using 2-stage implicit Radau-type RK
%
%   [t,Y] = MID_TERM_PC(B, y0, tspan, h) solves the linear ODE system
%   y' = B*y from tspan(1) to tspan(2) with time step h. Returns:
%       t : 1xN time vector
%       Y : solution matrix (m x N), columns = y at each time step
%
%   Optional argument:
%       doPlot : true/false (default = true) - plot solution if true
%
%   Example:
%       B = [-1 0 0; 0 -2 0; 0 0 -3];
%       y0 = [1;1;1];
%       tspan = [0 5];
%       h = 0.1;
%       [t,Y] = mid_term_pc(B, y0, tspan, h, true);

if nargin < 5
    doPlot = true;
end
if nargin < 4
    % Default test case
    B = [-1 0 0; 0 -2 0; 0 0 -3];
    y0 = [1;1;1];
    tspan = [0 5];
    h = 0.1;
end

% Extract time span
t0 = tspan(1);
tf = tspan(2);

% Time vector and steps
t = t0:h:tf;
N = length(t);

% Number of equations
m = length(y0);

% Preallocate solution matrix
Y = zeros(m, N);
Y(:,1) = y0;  % initial condition

% 2-stage Radau IIA coefficients
A = [5/12, -1/12;
     3/4,   1/4];
b = [3/4; 1/4];

% Precompute LU factorization for linear solve
M = build_radau_matrix(A, B, h);
[L,U,P] = lu(M);

% Time-stepping loop
for n = 1:N-1
    yn = Y(:,n);  % current solution

    % Build RHS for 2-stage RK
    F = [B*yn; B*yn];

    % Solve linear system M*K = F using LU
    K = U \ (L \ (P*F));

    % Extract stages
    k1 = K(1:m);
    k2 = K(m+1:end);

    % Update solution
    Y(:,n+1) = yn + h*(b(1)*k1 + b(2)*k2);
end

% Optional plotting
if doPlot
    figure
    plot(t, Y, 'LineWidth',1.5)
    xlabel('t')
    ylabel('y components')
    legend(arrayfun(@(i) sprintf('y%d',i),1:m,'UniformOutput',false))
    title('Radau 2-stage solution')
    grid on
end

end

%% Subfunction: build Radau RK matrix
function M = build_radau_matrix(A, B, h)
% BUILD_RADAU_MATRIX Constructs block matrix M for 2-stage Radau RK
%
%   M = BUILD_RADAU_MATRIX(A, B, h) returns the block matrix M such that
%   M*K = F represents the implicit RK system for stages K.
%
m = size(B,1);
I = eye(m);
M = [I - h*A(1,1)*B, -h*A(1,2)*B;
     -h*A(2,1)*B,    I - h*A(2,2)*B];
end