% Solutions to the heat flow equation 
% on one-dimensional rod length L
% Right end with Robin Condition
format compact;
L = 10;				% width of plate
Temp = 100;			% Constant temperature of rod, initially
tfin = 20;			% final time
k = 1;			    % heat coef of the medium
h = 0.5;            % Newton cooling constant
NptsX=151;          % number of x pts
NptsT=151;          % number of t pts
Nf=100;              % number of Fourier terms
x=linspace(0,L,NptsX);
t=linspace(0,tfin,NptsT);
[X,T]=meshgrid(x,t);

figure(1)
clf
a = zeros(1,Nf);
b = zeros(1,Nf);
U = zeros(NptsT,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)*exp(-k*(a(n))^2*T).*sin(a(n)*X);  % Temperature(n)
    U=U+Un;
    z0 = z0 + pi;
end
set(gca,'FontSize',[12]);
surf(X,T,U);
shading interp
colormap(jet)
xlabel('$x$','Fontsize',12,'interpreter','latex'); 
ylabel('$t$','Fontsize',12,'interpreter','latex'); 
zlabel('$u(x,t)$','Fontsize',12,'interpreter','latex');
axis tight
view([141 10])
print -depsc heat_rob1.eps		

figure(2)
clf

set(gca,'FontSize',[12]);
surf(X,T,U);
shading interp
colormap(jet)
view([0 90])             %create 2D color map of temperature
xlabel('$x$','Fontsize',12,'interpreter','latex'); 
ylabel('$t$','Fontsize',12,'interpreter','latex'); 
zlabel('$u(x,t)$','Fontsize',12,'interpreter','latex');
axis tight
colorbar
set(gca,'FontSize',[12]);
print -depsc heat_rob2.eps


