clear all

m=10;
a=zeros(m-1,m-1);
b=zeros(m-1,1);
T=zeros(m,1);
x=zeros(m,1);

k=14; h=100; W=0.02; t=0.005; Tb=373; T_inf=298;
L=0.1; 
delta_x=L/(m-1);

x(1)=0;
a(1,1)=-2-2*h*delta_x^2/k/t; a(1,2)=1;
b(1)=-2*h*delta_x^2/k/t*T_inf-Tb;
for i=2:m-2
    a(i,i-1)=1;
    a(i,i)=-2-2*h*delta_x^2/k/t;
    a(i,i+1)=1;
    b(i)=-2*h*delta_x^2/k/t*T_inf;
    x(i)=x(i-1)+delta_x;
end
a(m-1,m-2)=1; a(m-1,m-1)=-1;
b(m-1)=0;
x(m-1)=x(m-2)+delta_x; x(m)=L;

T(1)=Tb;
T(2:m,1)=inv(a)*b

% Exact solution
m=sqrt(h*2*W/(k*W*t));
x_e=0:0.001:L; s_x=size(x_e);
for i=1:s_x(2)
    T_e(i)=T_inf+(Tb-T_inf)*cosh(m*(L-x_e(i)))/cosh(m*L);
end

figure(1)
plot(x,T,'bv',x_e,T_e,'k-')
%plot(x,T,'ro')
grid on
xlabel('x (m)','fontname','times','fontangle','italic','fontsize',12)
ylabel('T (K)','fontname','times','fontangle','italic','fontsize',12)
set(gcf,'color',[1 1 1])
set(gca,'fontsize',12)
set(findobj('type','line'),'linewidth',2)
hold on


