Let me know if the looping approach I provided is helpful. By the way, your original proportional-only control law for xdot(2) could not drive the angular position to the desired θ reference of 1. So, I made a small modification to achieve that goal wonderfully. You can compare the code to study what changes were made.
kp = [1 2 4 8 16];
for j = 1:length(kp)
[t, x] = ode45(@(t, x) Time_domain_system(t, x, kp(j)), [0 20], [0; 0]);
plot(t, x(:,1)), hold on
end
hold off, grid, xlabel t
function xdot = Time_domain_system(t, x, kp)
theta_ref = 1;
m = 3;
L = 4;
J = 1/12*m*L^2;
c = 10;
g = 9.81;
k = 150;
% kp = 50;
% equivalent parameters
mx = J+m*L^2/4;
cx = c*L^2/4;
kx = k*L^2/4 - m*g*L/2;
xdot(1) = x(2);
% xdot(2) = kp*theta_ref/mx - (cx/mx)*x(2) - ((kx + kp )/mx)*x(1); % original
xdot(2) = kp*theta_ref - (cx/mx)*x(2) - ((kx + (mx*kp - kx))/mx)*x(1); % my proposal
xdot = xdot';
end
