%----------------------------------------------------------------------------
%  Script de calcul par differences finies 1D, stationnaire.
%  Resolution de l'equation de la chaleur :
%
%                      kd T,xx + f0 = 0
%
%  Conditions aux limites :
%                      T = T0 en x=0 et T = TL en x = L
%  
%  Forme discrete : alpha*T(i-1)+beta*T(i)+gamma*T(i+1)+fi=0 
%  
%----------------------------------------------------------------------------
clear all               % nettoyage de la memoire
close all               % fermeture des fenetres graphiques

%----- parametres physiques
L = 10;        		% longueur du domaine m
kd=2;		 	% coeff de conductivite w/°c-m
f0=15*0;               	% production     w/m3
T0=2;TL=100;		% conditions aux limites de Dirichlet

%----- parametres numeriques
nnt=input('Entrer le nombre de noeuds : ');
dx = L / (nnt - 1);     % pas de discretisation
k=zeros(nnt,nnt);       % initialisation de la matrice
f=zeros(nnt,1);         % initialisation du second membre
%
alpha=kd/dx^2; beta=-2*kd/dx^2; gamma=kd/dx^2; fi=f0;
%
for i=2:nnt-1           % schema de differences finies  [-1 2 -1]*kd/dx^2
     f(i) = -fi;
     k(i,i-1)=alpha;   
     k(i,i)=beta;
     k(i,i+1)=gamma;
end

%----- conditions aux limites de Dirichlet
k(1,1)=1;             f(1)=T0;   % en x=0
k(nnt,nnt)=1;         f(nnt)=TL; % en x=L

%----- resolution
T = k\f                 % calcul et affichage de la solution 

%----- affichage
x = 0:dx:L;             % vecteur des coordonnees des noeuds
hold off; plot(x,T','b -o')  % affichage

%------ solution analytique
d1=(TL-T0)/L + f0/(2*kd)*L;
d2=T0;
y = 0:0.1:L;
Texacte = (-.5*f0/kd)*y.^2+d1*y+d2;
hold on; plot(y,Texacte,'r');          % trace solution exacte
legend('Differences finies','Analytique');
xlabel('Abscisse x (m)');ylabel('Temperature (C)');
title('Distribution de la temperature dans la barre');