Perhaps you can try this approach.
G = tf(1, [1, 0, 0]);
opt = pidtuneOptions('DesignFocus', 'disturbance-rejection');
wc = 5;
C = pidtune(G, 'PID', wc, opt);
% Parameters
P.Xd = 3; % desired state for X(1)
P.Kp = C.Kp; % proportional gain
P.Ki = C.Ki; % integral gain
P.Kd = C.Kd; % derivative gain
% ODE function
function dX = f(t, X, P)
% error
e = P.Xd - X(1);
% PID controller
u = P.Kp*e + P.Ki*X(3) - P.Kd*X(2);
% disturbance
d = 1;
% original 2nd-order system
dX(1) = X(2);
dX(2) = d + u;
% X(3) is ∫ e dt
dX(3) = e;
% system dynamics
dX = [dX(1)
dX(2)
dX(3)];
end
% Solve and plot
[t, X] = ode45(@(t, X) f(t, X, P), [0, 10], zeros(3, 1));
plot(t, X(:,1)), grid on
title('Time response for state x_{1}')
xlabel('t')
ylabel('x_{1}')
