% Adveccion-Difusion 1D estacionaria por el Metodo de Volumenes Finitos
% Longitud del dominio
L=2;

% Difusividad
nu=1;

% Modulo de la velocidad
v=50/2.5;

% Condiciones de borde
phiI=0;
phiD=1;

% Cantidad de intervalos en los que sera discretizado el dominio
N=30;

% El valor de la seccion transversal es indistinto
S=1;

% Esquema de adveccion (UD, CD)
AS='CD';

% Calculo del paso de malla
h=L/N;

% Peclet de la malla
Pex=v*h/nu/2

% Calculo del volumen de los elementos
V=S*h;

% Ensamblado de la matriz difusiva
K=zeros(N,N);

for i=2:N-1
    K(i,i-1)=-nu*1/h;
    K(i,i)=nu*2/h;
    K(i,i+1)=-nu*1/h;
end

% Condiciones de borde
K(1,1)=3*nu/h;
K(1,2)=-nu/h;
K(N,N-1)=-nu/h;
K(N,N)=3*nu/h;

% Ensamblado de la matriz y el termino derecho para la adveccion
A=zeros(N,N);
b=zeros(N,1);

if (strcmp(AS,'UD'))
  for i=2:N-1
      A(i,i-1)=-v;
      A(i,i)=v;
  end
  % Condiciones de borde en la matriz
  A(1,1)=v;
  A(N,N-1)=-v;
  A(N,N)=v;
  % Condiciones de borde en el termino derecho
  b(1)=v*phiI;
elseif (strcmp(AS,'CD'))
  for i=2:N-1
    A(i,i-1)=-v/2;
    A(i,i+1)=v/2;
  end
  A(1,1)=v/2;
  A(1,2)=v/2;  
  A(N,N-1)=-v/2;
  A(N,N)=-v/2;
  b(1)=v*phiI;
  b(N)=-v*phiD;
else 
  disp('Esquema de adveccion no reconocido use UD ó CD')
end
  
% Incorporacion de las condiciones de borde de difusion
b(1)=b(1)+2*nu*phiI/h;
b(N)=b(N)+2*nu*phiD/h;


% Vector de abscisas
x=(h/2:h:(L-h/2))';
x=[0;x;L];

% Solucion
phi=(K+A)\b;

% Incorporacion de las condiciones de borde
phi=[phiI;phi;phiD];

% Graficacion 
hold off; %on
plot(x,phi,'r-');

% Graficacion de la solucion exacta
PeF=v*L/(2*nu)
phiEx=(exp(2*PeF*x/L)-1)/(exp(2*PeF)-1);

hold on;plot(x,phiEx,'b*-')

% Computo del error de la solucion numerica
e=norm(phi-phiEx)