% Conduction of heat in a semiinfinite slab with an initial temperature
% of zero and the x=0 boundary subject to a sudden change (step) of
% temperature from zero to T0. 
% 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
kappat=linspace(0,30,201);
kappat(1,1)=.001;
x=linspace(0,6,101);
Tdata=[];
for i=1:201
    T=ones(1,101)-erf(x./(2*sqrt(kappat(1,i))));
    Tdata=[Tdata;T];
end

figure(1);% Plot T(x,t) vs x for various t
xlabel('Distance, x')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
for k=1:5:201
    plot(x,Tdata(k,:),'k')
    hold on
    axis([0 6 -.2 1.2]);
    hold on
    %pause
end
text(.5,-.08,'\kappat = 0')
text(4 ,0.65,'\kappat = 30')
text(4,-.1,'\Delta(\kappat) = 0.6')
set(gca,'Box','on')
text(3,1.1,'Press Enter to Continue')
pause

figure(2)%Plot T(x,t) vs t for various x
for j=1:4:101
    plot(kappat, Tdata(:,j));
    hold on
    axis([0 30 -.2 1.2]);
    hold on
end
text(5,1.05, 'x = 0')
text(20, .3,'x = 6')
text(15,0,'\Deltax = 0.24')
xlabel('Scaled Time, \kappat')
ylabel('Dimensionless Temperature, T(x,t)/T_0')
text(15,1.1,'Press Enter to Continue')
pause

figure(3)%now do animation
axis([-1 6 -.2 1.2])
hold on
plot([0 0],[-.2,1.2]);
box on
xlabel('Distance, x')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
hold on
L=plot(x, Tdata(1,:),'k','EraseMode','xor');
hold on
texthandl=text(3, 1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=1:201
    set(L,'Ydata',Tdata(i,:))
    pause(.04)
end
hold on
set(texthandl,'String','Press Enter to Continue');
pause
    
figure(4)%3-D Plot of T(x,t) vs x and t
Tdata2=[];
for i=1:2:201
    Tdata2=[Tdata2;Tdata(i,:)];
end
axis([0 30 0 6 -.2 1.2])
hold on
kappat2=linspace(1,30,101);
[X,Y]=meshgrid(kappat2,x);
mesh(X,Y,Tdata2')
colormap([0 0 0])%cool)
text(25, 6, -.5,'Scaled Time,\kappat','Rotation',-30)
text(0, 4, -.5,'Distance, x','Rotation',13)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
grid on
text(30,5.7,1.1,'Press Enter to Continue')
view(240,30)