function S_N = FseriesOpt(f,T,t,N)

%Initialisera variabler
NoQpnt = 10*N; 			% Antal kvadraturpunkter
W = 2*pi/T;                     % Grundvinkelfrekvensen

dt = T/NoQpnt; 			% Delintervallens längd
qpnts = dt/2:T/NoQpnt:T-dt/2; 	% Räkna ut kvadraturpunkterna
ff = feval(f,qpnts); 		% Räkna ut f i kvaraturpunkter en
                                % gång för alla

% Axelgränser för plottning
ytop = max(abs(ff));
ymax = max(ff);
ymin = min(ff);

% Plotta f
subplot(121)
hold off
plot(t,feval(f,t),'r')
set(gca,'ylim',[ymin-1.25*ytop, ymax+1.25*ytop])
hold

S_N = ones(1,length(t))*mean(ff); % Räkna ut S_0=.5*a_0

% Plotta S_0:s graf och amplitud
subplot(121)
d=plot(t,S_N);
subplot(122)
hold off
plot([0,0],[0,mean(ff)])
set(gca,'xlim',[-1,N+1])
hold on
pause(.5)


for n=1:N
  a_n = 2/T*sum(ff.*cos(n*W*qpnts)*dt);      % Räkna ut a_n  
  b_n = 2/T*sum(ff.*sin(n*W*qpnts)*dt);      % Räkna ut b_n
  S_N = S_N + a_n*cos(n*W*t)+b_n*sin(n*W*t); % Addera till S_N
  subplot(121)
  delete(d);                         % Ta bort S_(N-1)
  d=plot(t,S_N);                     % Plotta S_N
  title(['n = ' num2str(n)])         % Plotta räknare
  subplot(122)
  plot([n,n],[0,sqrt(a_n^2+b_n^2)])  % Plotta amplituden
  pause(.5)
end
hold off





