%Conduction of heat in an infinite slab with initial constant 
%emperature T0 on (-L,L) and zero elsewhwere
% 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,3,91);
sn=[];
nterm=10;
ftime=1;
ntime=51;
tau=linspace(0,ftime,ntime);
tau(1,1)=.00001;
one=ones(1,91);
%T1=(ones(1,101)-xoverL);
framedata=[];
%framedata=[ones(1,31) zeros(1,60)];
for k=1:ntime
    den=2*sqrt(tau(1,k));
    V=0.5*(erf((one-xoverL)/den)+erf((one+xoverL)/den));
    framedata=[framedata;V];
end
framedata2=[fliplr(framedata(:,2:91)) framedata];
xoverL2=linspace(-3,3,181);

figure(1);% Plot T(x,t) vs x for various t

xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
for k=1:2:51
    plot(xoverL2,framedata2(k,:),'k')
    hold on
    axis([-3 3 -.1 1.2]);
    hold on
    %pause
end
%pause
text(-.95,.05,['\kappat/L^2 = ' num2str(0)])
%text(.16,.2,num2str(0.02))
text(-.2 ,.45,['\kappat/L^2 = ' num2str(1)])
%text(1.04,.8,'T = 0')
%text(-.15,.8,'T = T_0')
set(gca,'Box','on')
text(-2.5,1.1,'Press Enter to Continue')
pause

figure(2)%Plot T(x,t) vs t for various x
%plot(tau,ones(1,101))
hold on
for j=1:12:91
    plot(tau,framedata(:,j))
    hold on
    axis([0 1 0,1]);
    hold on
end
text(.4,.76,['x/L = ' num2str(0)])
text(.4,.62,['x/L = ' num2str(0.4)])
text(.4,.52,['x/L = ' num2str(.8)])
text(0.4,0.38,['x/L = ' num2str(1.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
axis([-3 3,-.1  1.2])
hold on
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
hold on
L=plot(xoverL2, framedata2(1,:),'k','EraseMode','xor');
hold on
texthandl=text(-2.5, 1.1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:51
    set(L,'Ydata',framedata2(i,:))
    pause(.05)
end
hold on
plot(xoverL2,framedata2(1,:),'k')
hold on
set(texthandl,'String','Press Enter to Continue');
pause


    
figure(5)%3-D Plot of T(x,t) vs x and t
[X,Y]=meshgrid(tau,(xoverL2));
mesh(X,Y,framedata2')
colormap(cool)
text(1,-1,-.4,'Distance, x/L','Rotation',12)
text(.1, -3,-.4,'Time, \kappat/L^2','Rotation',-32)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
view(60,30)
axis([0 1 -3 3  -.1 1.2])