Try this code, just an extension of your code.
clear variables
close all;
lim1 = -10;% no need to have t in pi
% units since omega is
% already measured in pi.
lim2 = 10;
M = 500; % length of
delta = (lim2-lim1)/(M-1);
t = lim1:delta:lim2;
A = 1.5;
f = 0.1;
omega= 2*pi*f;
phase= 0;
x = A*cos(omega*t+phase);
N = 10; % let say you want to average
% 10 terms of x
px = zeros(1,M);
for n = 1:N
px = px + (abs(A*cos(n*omega*t+phase))).^2;
end
px = (1/(2*N+1))*px;
%disp(px)
figure
plot(t,x,'linewidth',2);
hold on
plot(t,px,'linewidth',2)
a = ylabel('Signal');
set(a,'Fontsize',14);
a = xlabel('x');
set(a,'Fontsize',14);
a=title('x(t) and its average px');
legend('x(t)','px')
set(a,'Fontsize',16);
grid;
