clear
close all
L=1;                % Comprimento da haste (m)
T0=300;             % Temperatura na posição x=0 (K)
TL=1000;            % Temperatura na posição x=L (K)
TM=300;             % Temperatura para a equação da condição inicial (K)
a=0.1;              % Coeficiente de difusividade térmica (m^2/s)
nx=100;             % Número de intervalos na direção x
nt=20000;           % Número de intervalos no tempo
dx=L/nx;            % Tamanho dos elementos na direção x
dt=dx^2/(2*a); 	    % Intervalo de tempo (ajustado com a condição de estabilidade)
dt=0.25*dt;
T=zeros(nt,nx+1);   % Inicialização da matriz de temperatura
Ta=zeros(nt,nx+1);
erro=zeros(nt,nx+1);
t=zeros(nt+1,1);
for i=1:nx+1
    x(i)=(i-1)*dx;  % Cálculo das posições x ao longo da haste
    T(1,i)=TM*sin(pi/L*x(i))+(TL-T0)/L*x(i)+T0; % Condição inicial
    Ta(1,i)=TM*sin(pi/L*x(i))+(TL-T0)/L*x(i)+T0; % Solução analítica
    erro(1,i)=abs((T(1,i)-Ta(1,i))./Ta(1,i));
end
t(1)=0;
for n=1:nt
    n
    t(n+1)=n*dt;
    T(n+1,1)=T0;    % Temperatura no início da haste (condição de contorno)
    T(n+1,nx+1)=TL; % Temperatura no final da haste (condição de contorno)
    Ta(n+1,1)=T0;
    Ta(n+1,nx+1)=TL;
    for i=2:nx
        T(n+1,i)=T(n,i)+a*dt/(dx^2)*(T(n,i+1)-2*T(n,i)+T(n,i-1));   %Solução explícita
        Ta(n+1,i)=TM*exp(-a*(pi/L)^2*t(n+1))*sin(pi/L*x(i))+(TL-T0)/L*x(i)+T0; % Solução analítica
        erro(n+1,i)=abs((T(n+1,i)-Ta(n+1,i))./Ta(n+1,i));
    end
end
for n=1:100:nt+1
    n
    %plot(x,T(n,:))
    plot(x,100*erro(n,:))
    xlabel('x (m)')
    ylabel('erro em T (%)')
    %axis([0 L 0 max(max(T))])
    pause;
end
%pause
%plot(t,erro(:,40))