% Fourier series 
format compact;
L = 10;				% width of plate
Temp = 100;			% Constant temperature of rod, initially
h = 0.5;            % Newton cooling constant
NptsX=500;          % number of x pts
Nf=100;              % number of Fourier terms
X=linspace(0,L,NptsX);
a = zeros(1,Nf);
b = zeros(1,Nf);
U = zeros(1,NptsX);
U1 = zeros(1,NptsX);
U2 = zeros(1,NptsX);
U3 = zeros(1,NptsX);
z0 = 2.7;
for n=1:Nf
    z0 = fsolve(@(x) h*L*sin(x)+x*cos(x),z0);
    a(n) = z0/L;
    b(n)=(2*Temp*h/(a(n)*(L*h+(cos(a(n)*L))^2)))...
        *(1-cos(a(n)*L));  % Fourier coefficients
    Un = b(n)*sin(a(n)*X);  % Temperature(n)
    U = U+Un;
    if (n <= 5)
        U1 = U1+Un;
    end
    if (n <= 10)
        U2 = U2+Un;
    end
    if (n <= 20)
        U3 = U3+Un;
    end
    z0 = z0 + pi;
end
plot(X,U1,'m-','LineWidth',1.5);
hold on
plot(X,U2,'r-','LineWidth',1.5);
plot(X,U3,'-','Color',[0 0.5 0],'LineWidth',1.5);
plot(X,U,'b-','LineWidth',1.5);
plot([0 10],[100 100],'k-','LineWidth',1.5);
grid;
legend('n = 5','n = 10','n = 20','n = 100',...
    'location','southeast');
xlim([0 10]);
ylim([0 120]);
xlabel('$x$','Fontsize',12,'interpreter','latex'); 
ylabel('$f(x)$','Fontsize',12,'interpreter','latex');
set(gca,'FontSize',[12]);

print -depsc fourier_rob.eps		
