How can I plot the fourier series of a rectangular pulse by calculating the coefficients?

I want to find the fourier series of this recangular pulse shown in the code until N by using this formula

but I don't get the desired output.here is my code.could anyone tell what is the problem of my code?

here is the output I want to get

close all; clear; clc;
N = 2;
f = @(x) rectangularPulse(-1,1,x);
x = -2:0.01:2;
%2*p is the period
p = 1;
% the main function
plot(x,f(x),'LineWidth',2);
grid;
hold on;
grid minor;
xlim([-2 2]);
ylim([-0.1 1.1]);
x = linspace(-2,2,9).';
y = linspace(0,1,0.2);
a0 = (1/(2*p))*integral(f,-p,p);
an = zeros(1,N);
bn = zeros(1,N);
for n=1:N
    fan = @(x) rectangularPulse(-1,1,x).*cos((n*pi/p)*x);
    an(1,n) = (1/p)*integral(fan,-p,p);
    fbn = @(x) rectangularPulse(-1,1,x).*sin((n*pi/p)*x);
    bn(1,n) = (1/p)*integral(fbn,-p,p);
end
fs = a0 + sum(an.*cos((1:N).*x) + bn.*sin((1:N).*x),2);
plot(x,fs,'LineWidth',2);

Answer 1

[VBBV]

close all; clear; clc;
N = 2;
f = @(x) rectangularPulse(-1,1,x);
x = -2:0.001:2;
%2*p is the period
p = pi;
% the main function
plot(x,f(x),'LineWidth',2);
grid;
hold on;
grid minor;
xlim([-2 2]);
ylim([-0.1 1.1]);
x = linspace(-2,2,9).';
y = linspace(0,1,0.2);
a0 = (1/(2*p))*integral(f,-p,p);
an = zeros(1,N);
bn = zeros(1,N);
for n=1:N
    fan = @(x) rectangularPulse(-1,1,x).*cos((n*pi/p)*x);
    an(1,n) = (1/p)*integral(fan,-p,p);
    fbn = @(x) rectangularPulse(-1,1,x).*sin((n*pi/p)*x);
    bn(1,n) = (1/p)*integral(fbn,-p,p);
end
fs = a0 + sum(an.*cos((1:N).*x) + bn.*sin((1:N).*x),2);
plot(x,fs,'LineWidth',2);

[VBBV]

select a suitable period, for the series

[VBBV]

close all; clear; clc;
N = 2;
f = @(x) rectangularPulse(-1,1,x);
x = -2:0.001:2;
%2*p is the period
p = pi;
% the main function
plot(x,f(x),'LineWidth',2);
grid;
hold on;
grid minor;
xlim([-2 2]);
ylim([-0.1 1.1]);
x = linspace(-2,2,50).';
y = linspace(0,1,0.2);
a0 = (1/(2*p))*integral(f,-p,p);
an = zeros(1,N);
bn = zeros(1,N);
for n=1:N
    fan = @(x) rectangularPulse(-1,1,x).*cos((n*pi/p)*x);
    an(1,n) = (1/p)*integral(fan,-p,p);
    fbn = @(x) rectangularPulse(-1,1,x).*sin((n*pi/p)*x);
    bn(1,n) = (1/p)*integral(fbn,-p,p);
end
fs = a0 + sum(an.*cos((1:N).*x) + bn.*sin((1:N).*x),2);
plot(x,fs,'LineWidth',2);

[Niloufar]

thank you so much but if I change the value of N the graph becomes a little different for example for N = 6 the plot should be like this but my output is different.for N = 2 also has some differences

N = 6

[VBBV]

If you try with N = 10 or higher it looks similar and closely approximates the rectangular pulse

[VBBV]

Note that if you are comparing waveforms for N = 6 with other results as you shown, it may give different results. The output which you shown may have been obtained with a different upper limit summation for coefficients in the series.