[xheavenlyx]

Do one thing. When writing MATLAB code, divide each major segment into code blocks. As you have already used '%' for comment, use '%%' to create a block.

Run each block in part, one by one, when you reach the point where there is a problem, you can analyze the previous blocks output! Here is your code modified into code blocks:

%% Clean up
clear all
close all
clc

%% Define the system parameters; gain dimensions, time in minutes.
Kv = 5.0;
Tv = 10;
Tm = 0.2;

Kc = 1.3;
Tr = 5;
Td = 1;

L = 1; 

if Tr == 0
    Tr = 0.01;
end

%% Define the disturbances:time in minutes.
dr = 0; drt = 0; dd = -1; ddt = 0;
% Normal values of control signals: range 0-100%
r0 = 50; e0 = 0; u0 = 50; f0 = 50; d0 = 0; y0 = 50; m0 = 50;
 
% Initialise simulation variables.
f = f0; ia = 0; m = m0; mp = m0; t = 0; u = u0; y = y0;

%
% Establish step length an run time.
dt = Tv/200;
stop = 800;

%
% Initialise vectors for display purposes.
T(1) = 0;
Y(1) = y0;

%% Establish a process delay.
n = round(L/dt);
for i = 1:n
    Q(i) = f0;
end

%%
% Enter control parameters.
echo off;
disp (['classical form of PID with derivative feedback']);

% Controller gains.
Ki = Kc*dt/Tr; 
Kd = Kc*Td/dt;

%%
% Recursion of loop dynamics.
for j = 1:stop
    t = t+dt;
    T(j+1) = t;
    
    %
    % Dynamics of controller.
    if t < drt
        r = r0;
    else
        r = r0+dr;
    end
    R(j+1) = r;
    e = r-m;
    % Apply anti-windup to I action.
    if u == 0
        u = u0+Kc*e+ia-Kd*(m-mp);
    elseif u == 100
        u = 0+Kc*e+ia-Kd*(m-mp);
    else
        ia = ia+Ki*e;
        u = u0+Kc*e+ia-Kd*(m-mp);
    end
    % Force controller Saturation.
    if u < 0
        u = 0;
    end
    if u > 100
        u = 100;
    end
    U(j+1) = u;
    
%%
    % Dynamics of i/p converter, actuator, positioner and valve.
    f = f+(Kv*(u-u0)-(f-f0))/Tv*dt;

    % Force process saturation.
    if f < 0 
        f = 0;
    end
    if f > 100
        f = 100;
    end
    
    %
    % Dynamics of delay.
    if t < ddt
        d = d0;
    else
        d = d0+dd;
    end
    y = Q(1);
    for i = 1:n-1
        Q(i) = Q(i+1);
        Q(n) = f+d;
    end
    Y(j+1) = y;
    
    %
    % Dynamics of sensor, transducer and transmitter
    mp=m;
    m = m+(1*(y-y0)-(m-m0))/Tm*dt;
    %
    
end
%%
plot (T,Y,'b')
grid on
xlabel('Time')
ylabel('Simulation output')
title('Simulation Results','FontWeight','bold')
legend('System Output','location','SouthEast') 
zoom
%

For example here after you complete running block with "Establish a process delay" : Int he command window type "Q" to get its values:

>> Q

Q =

Columns 1 through 13

50 50 50 50 50 50 50 50 50 50 50 50 50

Columns 14 through 20

50 50 50 50 50 50 50

>>

And you know, I could not get exactly what the problem is, cuz I was getting the graph. The Lc/Tc I cannot understand. I had Chemical process control the last sem as an elective. But we did not go in depth Please explain the process a bit in your next post if possible with the results of this method (code blocks).