Plotting a system of fractional order differential equations

Hi. I have been trying to plot fractional order differential equations with MATLAB. I have tried to use this technique called "Predictor-Corrector Method" and got a code listing using AI. The code seems to be giving the plots, though they do not match with the expected results. I will attach the code listing below. But my question is how to actually plot such systems? Is there a standard method in-built, or a widely accepted add-on/package? Thanks. ```% Parameters alpha = 0.9; LAMBDA = 3; mu = 0.4; beta = 0.004; gamma = 0.0005; m = 0.03; n = 0.006; delta = 0.02; rho = 0.003; h = 0.1; % Time step T = 300; % Total time

% Time steps time = 0:h:T; time_steps = length(time);

% Initial conditions s = zeros(1, time_steps); i = zeros(1, time_steps); r = zeros(1, time_steps); s(1) = 85; i(1) = 3; r(1) = 0;

% Predictor-Corrector Method for k = 1:(time_steps - 1) % Predictor step (Adams-Bashforth) if k == 1 ds_pred = h * (LAMBDA - mu * s(k) - (beta * s(k) * i(k)) / (1 + gamma * i(k))); di_pred = h * ((beta * s(k) * i(k)) / (1 + gamma * i(k)) - (m * i(k)) / (1 + n * i(k)) - (mu + delta + rho) * i(k)); dr_pred = h * ((m * i(k)) / (1 + n * i(k)) + delta * i(k) - mu * r(k)); else ds_pred = h / 2 * ((LAMBDA - mu * s(k) - (beta * s(k) * i(k)) / (1 + gamma * i(k))) + ... (LAMBDA - mu * s(k-1) - (beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1)))); di_pred = h / 2 * (((beta * s(k) * i(k)) / (1 + gamma * i(k)) - (m * i(k)) / (1 + n * i(k)) - (mu + delta + rho) * i(k)) + ... ((beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1)) - (m * i(k-1)) / (1 + n * i(k-1)) - (mu + delta + rho) * i(k-1))); dr_pred = h / 2 * (((m * i(k)) / (1 + n * i(k)) + delta * i(k) - mu * r(k)) + ... ((m * i(k-1)) / (1 + n * i(k-1)) + delta * i(k-1) - mu * r(k-1))); end

    % Update predicted values
    s(k + 1) = s(k) + ds_pred;
    i(k + 1) = i(k) + di_pred;
    r(k + 1) = r(k) + dr_pred;

    % Corrector step (Adams-Moulton)
    ds_corr = h / 2 * ((LAMBDA - mu * s(k + 1) - (beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1))) + ...
                       (LAMBDA - mu * s(k) - (beta * s(k) * i(k)) / (1 + gamma * i(k))));
    di_corr = h / 2 * (((beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1)) - (m * i(k + 1)) / (1 + n * i(k + 1)) - (mu + delta + rho) * i(k + 1)) + ...
                       ((beta * s(k) * i(k)) / (1 + gamma * i(k)) - (m * i(k)) / (1 + n * i(k)) - (mu + delta + rho) * i(k)));
    dr_corr = h / 2 * (((m * i(k + 1)) / (1 + n * i(k + 1)) + delta * i(k + 1) - mu * r(k + 1)) + ...
                       ((m * i(k)) / (1 + n * i(k)) + delta * i(k) - mu * r(k)));

    % Update corrected values
    s(k + 1) = s(k) + ds_corr;
    i(k + 1) = i(k) + di_corr;
    r(k + 1) = r(k) + dr_corr;
end

% Plot the results figure; plot(time, s, '-o', 'DisplayName', 's(t)'); hold on; plot(time, i, '-s', 'DisplayName', 'i(t)'); plot(time, r, '-d', 'DisplayName', 'r(t)'); xlabel('Time'); ylabel('Population'); title('Solution of the System of Fractional Differential Equations'); legend; grid on;

% Set plot range xlim([0 30]); ylim([0 100]);```