% Conduction of heat in a slab with left face at T0 and right face at zero
% and slab initially at zero temperature. The problem is solved 
% approximately by finite differences in space and state variables in time.
% The diffusivity is kappa=k/(rho*c)
% 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);
figure(1);clf
figure(2);clf
figure(3);clf
figure(4);clf
figure(5);clf
xoverL=linspace(0,1,21);
N=19;% 19 internal nodes
A=zeros(19,19);
B=zeros(19,1);
B(1,1)=400;
for i=1:19
    A(i,i)=-800;
end
for i=1:18
    A(i,i+1)=400;
    A(i+1,i)=400;
end

T=.001;
[Ad,Bd]=c2d(A,B,T);
x=zeros(19,1);
framedata=x';
for k=1:200
    x=Ad*x+Bd;
    framedata=[framedata;x'];
end
framedata=[ones(201,1) framedata zeros(201,1)];
framedata(1,1)=0;


figure(1);% Plot T(x,t) vs x for various t
f1=[ 0 0.08 .2 .31 .43 .51 .58 .64 .76 .88 1];
xpatch1=[f1 ones(1,11)-f1 f1(1,1)];
y1=[-.1 -.04 -.16 -.12 -.05 -.14 -.07 -.11 -.16 -.12 -.08];
ypatch1=[y1 1.1*ones(1,11)+(y1) y1(1,1)];
patch(xpatch1,ypatch1,'c')
hold on
plot(xoverL,framedata(1,:),'k')
hold on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
for k=11:20:201
    plot(xoverL,framedata(k,:),'k')
    hold on
    axis([-.2 1.2 -.2 1.2])
end
text(.05,.05,['\kappat/L^2 = ' num2str(0)])
text(.16,.2,num2str(0.02))
text(.4,.54,num2str(0.2))
text(1.04,.8,'T = 0')
text(-.15,.8,'T = T_0')
set(gca,'Box','on')
text(.1,1.13,'Press Enter to Continue')
pause

figure(2)%Plot T(x,t) vs t for various x
tau=linspace(0,.2,201);
plot(tau,ones(1,201))
hold on
for j=3:2:21
    plot(tau,framedata(:,j)')
    hold on
end
text(.02,.97,['x/L = ' num2str(0)])
text(.03,.8,['x/L = ' num2str(0.1)])
text(.1,-.08,['x/L = ' num2str(1)])
text(0.04,0.6,['x/L = ' num2str(.2)])
xlabel('Dimensionless Time, \kappat/L^2')
ylabel('Dimensionless Temperature, T(x,t)/T_0')
text(.02,1.1,'Press Enter to Continue')
pause

figure(3)%now do animation
xpatch=[0 0 .2 .3 .5 .6 .85 1 1 .82 .66 .43 .2 0];
ypatch=[-.1 1.1 1.15 1.05 1.13 1.1 1.14 1.1 -.1 -.03 -.14 -.1 -.16 -.1];
patch(xpatch,ypatch,'c')
axis([-.2 1.2 -.2 1.2])
hold on
text(-.15,.8,'T = T_0')
text(1.05,.8,'T = 0')
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
hold on
L=plot(xoverL, framedata(1,:),'k','EraseMode','xor');
hold on
texthandl=text(.4,1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:201
    set(L,'Ydata',framedata(i,:))
    pause(.02)
end
hold on
plot(xoverL,framedata(1,:),'k')
hold on
set(texthandl,'String','Press Enter to Continue');
pause

figure(4)%Color plot
axis([-.2 1.2 -.2 1.2 ]);
hold on
framedata(:,1)=ones(201,1);
colormap(cool);
zdata=zeros(21,21);
for j=1:21
    zdata(j,:)=framedata(1,:);
end
Q=pcolor(xoverL,xoverL,zdata);
xlabel('Dimensionless Distance, x/L');
hold on
ylabel('Dimension y');
texthandl2=text(0,1.1,'Press Enter to Continue');
pause

set(texthandl2,'String','        ');
for i=2:2:200
    for j=1:21
        zdata(j,:)=framedata(i,:);
    end
    Q=pcolor(xoverL,xoverL,zdata);
    pause(.005)
end
set(texthandl2,'String','Press Enter to Continue');
pause
    
figure(5)%3-D Plot of T(x,t) vs x and t
[X,Y]=meshgrid(xoverL,tau);
mesh(X,Y,framedata)
colormap(cool)
xlabel('Dimensionless Distance, x/L','Rotation',-34)
ylabel('Dimensionless Time, \kappat/L^2','Rotation',12)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
view(60,30)
axis([0 1 0 .2 0 1])