% script to plot Fourier series representation of triangle function
% M is number of Fourier terms
% N is number of x values in [0,1]
% M and N must be supplied from outside
alpha=0.5;
x=linspace(0,1,N); %x is the independent variable
ysteady=(1+1/2/alpha/alpha)*x-x.*x/2/alpha/alpha;

m=1:M;  %m labels the Fourier terms
FourierCoef=(2/pi./m+4/pi/pi/pi/alpha/alpha./m./m./m).*mod(m,2)...
   -2/pi./m.*(~mod(m,2)); %Fourier decomposition of ysteady
SineFunctions=sin(pi*m'*x);  %each row of this matrix is one of the sine functions
ytrans=-(FourierCoef.*exp(-m.*m*alpha*alpha*pi*pi*t))*SineFunctions; %transient part
y=ysteady+ytrans; %the complete solution
plot(x,y)