How to illustrate the dynamic system output for a PID controller system?

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 when I try to illustrate the x3 output, it just scrambles everything and after 26th iteration it doesn't even evaluate(bc numbers too big or small).It should control my output in a specified range with my reference, control and error signal (r(n), u(n), e(n) respectively) but i think it doesn' do that. I'm open to any help or idea. Thanks in advance.

clear; clc; close all;

% Initialization
Kp0 = 1.5;

Ki0 = 0.1;

Kd0 = 0.25;

u=0:1;

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

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

d2 = 1;

t = 1:20;
r = 0.45 + 0.15*sin(2*pi*(1/200)*t);

for n = 3:length(r)
    % Compute error
    e(n) = r(n) - x3(n);
     
    % Update PID gains
    Kp(n+1) = Kp0 + e(n);
    Ki(n+1) = Ki0 + e(n);
    Kd(n+1) = Kd0 + e(n);

    % 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) > max(u)
        u(n) = max(u);
    elseif u(n) < min(u)
        u(n) = min(u);
    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);
   

    % Plotting
    figure(1);

    subplot(3,1,1);
    plot(x3, 'r');
    title('x3 - r');

    subplot(3,1,2);
    plot(r);

end

Pseudocode

What I need to do

I try to illustrate the x3 output, it just scrambles everything and after 26th iteration it doesn't even evaluate(bc numbers too big or small).