The easiest way would be to use fplot(). Also, use assume() instead of evalin
syms x k L n
% evalin(symengine,'assume(k,Type::Integer)');
assume(k,'integer')
assume(n,'integer')
%Obtain an and bn
a = @(f,x,k,L) int(f*cos(k*pi*x/L)/L,x,-L,L);
b = @(f,x,k,L) int(f*sin(k*pi*x/L)/L,x,-L,L);
%Fourier Series is calculated
fs = @(f,x,n,L) a(f,x,0,L)/2 + ... %a0/2 a0=an in the interation number 0 [k or n = 0]
symsum(a(f,x,k,L)*cos(k*pi*x/L) + b(f,x,k,L)*sin(k*pi*x/L),k,1,n); %+an*cos(nwt)+bn*sin(nwt)
f = x %Function
pretty(fs(f,x,5,pi)); %fs(f(t), x, n, T) ----- fs(function, variable x, number of terms, period)
sys_sum = fs(f,x,5,pi);
fplot(sys_sum)
