% This module solves the conduction in a finite slab with convective
% boundaries.  The temperature at the left boundary layer is suddenly elevated
% T0 degrees above the original equilibrium temperature.  The Biot Number at
% both boundaries is assumed the same (unity). 
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
%Bi=.1;
%lambdaL=[.44352 3.203994 6.314854 9.445951 12.58227 15.72069 18.86016];
%lambdaL=[lambdaL 22.00024 25.1407 28.2814 31.42229 34.56331 37.70442];
Bi=1;
lambdaL=[1.306542 3.673195 6.58462 9.631684 12.72324 15.83411];
lambdaL=[lambdaL 18.95497 22.08166 25.21203 28.34486 31.47944 34.61528];
xoverL=linspace(0,1,201);
delta=.005;%delta(x/L), integration increment
for i=1:12
    z=lambdaL(1,i)*xoverL;
    phi(i,:)=(Bi/lambdaL(1,i))*sin(z)+cos(z);
end
%Check Orthogonality of Eigenfunctions
for i=1:12
    for j=1:12
        I(i,j)=(sum(phi(i,:).*phi(j,:))-(phi(i,1)*phi(j,1)+phi(i,201)*phi(j,201))/2)*delta;
    end
end
den=Bi+2;
A=(Bi+1)/den;
B=-Bi/den;
f=A*ones(1,201)+B*(xoverL);% Steady State Temperature in Slab
sum1=zeros(1,201);
for i=1:12
    b(1,i)=(1/I(i,i))*(sum(-f.*phi(i,:))-(-f(1,1)*phi(i,1)-f(1,201)*phi(i,201))/2)*delta;
    sum1=sum1+b(1,i)*phi(i,:);
end
T=sum1+f;
tempdata=[T];
T(1,1)=0;
tau=linspace(0,2,201);% tau=kappa*t/(L^2), Dimensionless Time
for t=2:201
    sum2=zeros(1,201);
    for i=1:12
        sum2=sum2+b(1,i)*exp(-(lambdaL(1,i)^2)*tau(1,t))*phi(i,:);
    end
    T=sum2+f;
    tempdata=[tempdata;T];
end

figure(1);clf% Plot T vs time for 11 values of x/L
axis([0 2 -.1 1.1]);
hold on
for i=1:20:201
    plot(tau,tempdata(:,i));
    hold on
end
box on
ylabel('Dimensionless Temperature Change, T(x,t)/T_0')
xlabel('Dimensionless Time, \kappat/L^2')
text(1.7,.28,'x/L = 0');
text(1.7,.7,'x/L = 1');
text(1.7,.1,'\Deltax/L = 0.1');
text(1, 1,'Press Enter to Continue');
pause
xoverL1=[-.1 -.1 xoverL 1.1];
tempdata=[ones(201,1) ones(201,1) tempdata zeros(201,1)];
tempdata(1,:)=zeros(1,204);
tempdata(1,1)=1;

figure(2);clf% Plot T versus x/L for various times
patch([0 0 1 1 0],[-.1 1.1 1.1 -.1 -.1],'c')
hold on
axis([-.2 1.2 -.1 1.1])
hold on
for t=1:10:101
    plot(xoverL1,tempdata(t,:))
    hold on
end
text(.6,.9,['Biot No. = ' num2str(Bi)])
xlabel('Dimensionless distance, x/L')
ylabel('Dimensionless Temperature Change, T(x,t)/T_0')
text(.2,.025,'\kappat/L^2 = 0')
text(.2,.1,'0.1')
text(.2,.53,'1')
box on
text(.1, 1,'Press Enter to Continue')
pause

figure(3);clf% Animate T vs x/L for closely spaced values of time
axis([-.2 1.2 -.1 1.1])
hold on
patch([0 0 1 1 0],[-.1 1.1 1.1 -.1 -.1],'c')
hold on
L=plot(xoverL1,tempdata(1,:),'EraseMode','xor');
ylabel('Dimensionless Temperature Change, T(x,t)/T_0')
xlabel('Dimensionless Distance, x/L')
text(.6,.9,['Biot No. = ' num2str(Bi)])
plot([-.2 1.2],[-.005 -.005])
hold on
box on
T1=text(.1,1,'Press Enter to Animate');
pause
set(T1,'String','                     ')
for t=1:201
    set(L,'Ydata',tempdata(t,:));
    pause(.1)
end
box on
set(T1,'String','Press Enter to Continue');
pause
figure(4);clf
[X,Y]=meshgrid(xoverL1,tau);
mesh(X,Y,tempdata);
colormap(cool);
axis([-.1 1.1 0 2 0 1.1]);
hold on
text(0,0,-.2,'Location, x/L','Rotation',-34)
text(1.1,.4,-.2,'Dimensionless Time, \kappat/L^2','Rotation',12)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
view(60,30)
text(0.1,2,.7,'Press Enter to Continue');
pause