can this code be improved?

I would like to know if this code is well written and if it can be improved? any help would be appreciated. Thanks.

% --- Problem Setup ---
n_e = 100;                % 100 ODEs
t_span = [0, 0.2];      % Reduced time for stability demo
h = 1e-4;               % Step size 
t_steps = 0:h:t_span(2);
dx = 2/99;              
xx = 0:dx:2;            
y0 = ones(n_e, 1);        % Initial condition y(0) = [1, ..., 1]^T


%% 1. Radau IIA (Implicit Collocation - s=2)
% Butcher Tableau for Radau IIA (Order 3)

A_rad = [5/12, -1/12; 3/4, 1/4];
b_rad = [3/4, 1/4];
y_rad = y0;
imp_rad = zeros(n_e, length(t_steps));
imp_rad(:,1) = y0;

fprintf('Running Radau IIA...\n');

Running Radau IIA...

for k = 1:length(t_steps)-1
    y_old = imp_rad(:,k);
    options = optimoptions('fsolve','Display','off');

    % Solve Stage Equations (Zi = y_old + h*sum(a_ij*f(Zj)))
    Z_guess = [y_old; y_old];
    Z_sol = fsolve(@(Z) radau_stages(Z, y_old, h, A_rad, n_e), Z_guess, options);
    
    Z1 = Z_sol(1:n_e); Z2 = Z_sol(n_e+1:end);
    b_rad = [3/4, 1/4];

    imp_rad(:,k+1) = y_old + h*( ...
    b_rad(1)*combustion_fun(Z1) + ...
    b_rad(2)*combustion_fun(Z2) );
end

%% 2. Adams-Bashforth
% Formula: y_n = y_{n-1} + h * (1.5*f_{n-1} - 0.5*f_{n-2})

y_ab = y0;
imp_ab = zeros(n_e, length(t_steps));
imp_ab(:,1) = y0;

fprintf('Running Adams-Bashforth 2...\n');

Running Adams-Bashforth 2...

% Step 1: Forward Euler 
f_minus_1 = combustion_fun(y0);
imp_ab(:,2) = y0 + h * f_minus_1;

for k = 2:length(t_steps)-1
    f_curr = combustion_fun(imp_ab(:,k));

    % AB Formula
    imp_ab(:,k+1) = imp_ab(:,k) + h * (1.5 * f_curr - 0.5 * f_minus_1);
    f_minus_1 = f_curr;
    
    if any(isnan(imp_ab(:,k+1))), error('AB2 unstable. Reduce h.'); end
end

figure('Position', [100, 100, 800, 400]); 
subplot(1,2,1); hold on;
plot_indices = round(linspace(1, length(t_steps), 6));
for idx = plot_indices
    plot(xx, imp_rad(:,idx), 'DisplayName', sprintf('t=%.2f', t_steps(idx)));
end
title('Radau IIA (Implicit - Stable)'); xlabel('xx'); ylabel('y(t)'); legend; grid on;

subplot(1,2,2); hold on;
for idx = plot_indices
    plot(xx, imp_ab(:,idx), '--', 'DisplayName', sprintf('t=%.2f', t_steps(idx)));
end
title('Adams-Bashforth 2 (Explicit)'); xlabel('xx'); ylabel('y(t)'); legend; grid on;

%sgtitle('Transition of Vector y(t) (Test 1)');

%% --- Functions ---
function F = combustion_fun(y) % Source: [2]
    n = length(y); F = zeros(n, 1);
    R=5; delta=20; alpha=1; D=R*exp(delta)/delta; T0=1+alpha;
    dx_sq = (2/(n-1))^2;
    for i = 2:n-1
        F(i) = (y(i+1)-2*y(i)+y(i-1))/dx_sq + D*(T0-y(i))*exp(-delta/y(i));
    end
    F(1) = 0; F(n) = 0; % Boundary Conditions [1]
end

function Res = radau_stages(Z, y_old, h, A, n)
    Z1 = Z(1:n); Z2 = Z(n+1:end);
    f1 = combustion_fun(Z1); f2 = combustion_fun(Z2);
    Res = [Z1 - y_old - h*(A(1,1)*f1 + A(1,2)*f2);
           Z2 - y_old - h*(A(2,1)*f1 + A(2,2)*f2)];
end

Answer 1

[John D'Errico]

You have asked several questions. 1. Is it well written? 2. Can it be improved? Something you did not ask is if it SHOULD be improved. That is, is it worth wasting time to improve code that works, and at least, apparently does what you need to do. Does it run in a reasonable time? If so, move onto the next thing you need to do. Eventually, your coding skills will improve, and you will be writing better code anyway.

To your questions as posed...

Is it well written? It is not terrible code. Not even bad, certainly compared to much of what we see online. You have comments, making it possible to sort of read the code. That is good. You use variable names that sort of make sense, making your code readable. Again, good. Overall, on a scale of A to F, with A being the best grade, I'd give the code a B. It could depend on how other students were answering the same question, and admittedly, many students write complete dreck for code these days. Also, I've not checked to see if the code works. If it did not run correctly, that would be a significant downside.

I'd not call it good code though. You have your functions hard coded inside the solvers. That is a poor coding style, since tomorrow you will want to solve a different problem. You have boundary conditions hard coded into the solver, which is again a poor programming paradigm. Learn to separate solvers from the current problem. Pass in an objective into a solver.

Better yet, learn to use existing solvers to solve your problems. Since this is probably a homework assignment, that is not terribly pertinent, but it is something you need to learn. You should never be writing your own ODE solver code, when high quality tools already exist for your use. But since your long term goal will be to learn to use the existing solvers in MATLAB, you need to learn to write functions, and to learn to pass functions to other functions. This is why I do not really like your code, at least why it is not A level code in my eyes.

One thing you could consider doing is to look at tools like the ODE solvers in MATLAB. Think about their basic interface, where you pass in the right hand side for a differential equation or a system of equations. This is basically the same idea embodied in fsolve, which you do know how to use.

Can your code be improved? Well, yes. I just pointed out several ways you might have done so. As I also said, I see no reason to spend the time to perfect working, viable code. Go on to the next problem you need to solve.

Answer 2

[Torsten]

Better use separate functions for problem setup, ODE solution and plotting and call these three functions from one main driver. Especially for larger problems than yours from above, this will help to intuitively understand the structure of the code.

In previous questions, you coded a stepsize control. Don't you have to use it in the code from above ? Especially if "fsolve" for RADAU II a doesn't converge, you will have to reduce the stepsize and try again.

And are you sure you can use MATLAB's "fsolve" for the arising nonlinear systems ? If you try to code RADAU II a on your own, usually writing a solver for the nonlinear equations would also be a task of yours.

One little thing at the end: If "imp" means "implicit", I don't understand why you call the solution obtained from Adams-Bashford "imp_ab". As you state yourself, Adams-Bashford is explicit.