clear all;
close all;

x = linspace(-3,6,1000);

b0 = 1;
soln = b0*ones(1,length(x));

for n=1:4
   an = (1-2*cos(4*n*pi/3)+cos(2*n*pi/3))/(n*pi)
   bn = (2*sin(4*n*pi/3)-sin(2*n*pi/3))/(n*pi)
   soln = soln + an*sin(2*n*pi*x/3) + bn*cos(2*n*pi*x/3);
end

plot(x,soln);
grid on;
hold on;
for n=5:10
   an = (1-2*cos(4*n*pi/3)+cos(2*n*pi/3))/(n*pi)
   bn = (2*sin(4*n*pi/3)-sin(2*n*pi/3))/(n*pi)
   soln = soln + an*sin(2*n*pi*x/3) + bn*cos(2*n*pi*x/3);
end
plot(x,soln,'r');

for n=11:15
   an = (1-2*cos(4*n*pi/3)+cos(2*n*pi/3))/(n*pi)
   bn = (2*sin(4*n*pi/3)-sin(2*n*pi/3))/(n*pi)
   soln = soln + an*sin(2*n*pi*x/3) + bn*cos(2*n*pi*x/3);
end
plot(x,soln,'m');

for n=16:75
   an = (1-2*cos(4*n*pi/3)+cos(2*n*pi/3))/(n*pi)
   bn = (2*sin(4*n*pi/3)-sin(2*n*pi/3))/(n*pi)
   soln = soln + an*sin(2*n*pi*x/3) + bn*cos(2*n*pi*x/3);
end
figure;
plot(x,soln,'k');


