How to properly iterate functions/equations

Question

I'm trying to make a PID Controller with cstr in it for my homework in matlab environment, I've managed to come up with something promising but encountered an index error. I appreciate any ideas, thanks in advance.

Pseudocode, System Dynamics and what I need to do next

clear; clc; close all;

% Initialization
Kp0 = 0.75;
Kp(3) = Kp0;

Ki0 = 0.1;
Ki(3) = Ki0;

Kd0 = 0.25;
Kd(3) = Kd0;

umax = 1;
umin = 0;
u(2) = 0.1;

% CSTR system parameters
Ts = 0.1;
Da1 = 3;
Da2 = 0.5;
Da3 = 1;

x1(1) = 0.5;
x2(1) = 0.5;
x3(1) = 1.5;

d2 = 1;

for n = 3:100
    % Compute error
    e(n-2) = x1(n-2) - x3(n-2);
    
    % PID Control Signal
    u(n) = u(n-1) + Kp(n) * (e(n) - e(n-1)) + Ki(n) * e(n) + Kd(n) * (e(n) - 2*e(n-1) + e(n-2));
    
    % Saturate Control Signal
    if u(n) > umax
        u(n) = umax;
    elseif u(n) < umin
        u(n) = umin;
    end
    
    % Apply u to CSTR system
    k1x1 = Ts*f1_func(x1(n), x2(n));
    k1x2 = Ts*f2_func(x1(n), x2(n), d2, u(n));
    k1x3 = Ts*f3_func(x2(n), x3(n), d2);

    k2x1 = Ts*f1_func(x1(n) + k1x1/2, x2(n) + k1x2/2);
    k2x2 = Ts*f2_func(x1(n) + k1x1/2, x2(n) + k1x2/2, d2, u(n));
    k2x3 = Ts*f3_func(x2(n) + k1x2/2, x3(n) + k1x3/2, d2);

    k3x1 = Ts*f1_func(x1(n) + k2x1/2, x2(n) + k2x2/2);
    k3x2 = Ts*f2_func(x1(n) + k2x1/2, x2(n) + k2x2/2, d2, u(n));
    k3x3 = Ts*f3_func(x2(n) + k2x2/2, x3(n) + k2x3/2, d2);

    k4x1 = Ts*f1_func(x1(n) + k3x1, x2(n) + k3x2);
    k4x2 = Ts*f2_func(x1(n) + k3x1, x2(n) + k3x2, d2, u(n));
    k4x3 = Ts*f3_func(x2(n) + k3x2, x3(n) + k3x3, d2);

    x1(n + 1) = x1(n) + 1/6 * (k1x1 + 2*k2x1 + 2*k3x1 + k4x1);
    x2(n + 1) = x2(n) + 1/6 * (k1x2 + 2*k2x2 + 2*k3x2 + k4x2);
    x3(n + 1) = x3(n) + 1/6 * (k1x3 + 2*k2x3 + 2*k3x3 + k4x3);

    % Update PID gains
    Kp(n+1) = Kp0 + e(n);
    Ki(n+1) = Ki0 + e(n);
    Kd(n+1) = Kd0 + e(n);

    % Plotting
    figure(1);

    subplot(3,1,1);
    plot(n, x1(n), '-o');
    title('x1');
    hold on;

    subplot(3,1,2);
    plot(n, x2(n), '-o');
    title('x2');
    hold on;

    subplot(3,1,3);
    plot(n, x3(n), '-o');
    title('x3');
    hold on;
end

I've tried to implement given pseudocode just as it is. but with given functions, my script is trying to index an element that has not yet been created in my equation for u(n). (e(n))