% difusionrobindf1dimpl.m
%
% Programita para resolver la ecuacion de Difusion
% con condiciones de Dirichlet en x=0 y Robin en x=L.
% Metodo: Diferencias finitas en espacio y Euler IMPLICITO en tiempo.
%
%  u_t(x,t) - k u_xx(x,t) = f(x,t),   0 < x < L,  t > 0,
%  u(0,t) = a(t),    k u'(L,t) + H1 u(L,t) = H2 uE(t),  t > 0,
%  u(x,0) = u0(x),   0 < x < L.

%% Parametros del problema
L = 1;  T = 2.5;  k = 1;
f = @(x,t)(10*exp(-100*(x-.5).^2)*(t<1.5));
a = @(t)(0);
H1 = 3;  H2 = 3;  uE = @(t)(0);
u0 = @(x)(zeros(size(x)));

%% Parametros del metodo de resolucion
N = 30;
h = L/N;
lambda = 0.25; % la mitad de lo necesario para tener estabilidad
deltat = lambda*h^2/k; 

%% Condicion inicial
X = linspace(0,L,N+1)';
U = u0(X);
% le damos valor a la U en el nodo ficticio
U(N+2) = (H2*uE(0)-H1*U(N+1))*2*h/k + U(N);

ejes = [0 L 0 1];
figure(1);  plot(X,U(1:N+1)); 
title('t = 0');
axis(ejes);  pause

%% Armado de la matriz
unos = ones(N+2,1);
columnas = [-lambda*unos (1+2*lambda)*unos -lambda*unos];

matriz = spdiags(columnas, [-1 0 1], N+2, N+2);
matriz(1,1:3) = [1 0 0];
matriz(N+2, N:N+2) = [-1 2*h*H1/k 1];

%% Resolucion
F=zeros(N+2,1);

for t=deltat:deltat:T
    % Ensamblado del lado derecho
    U(1) = a(t);
    U(N+2) = 2*H2*h/k*uE(t);
    
    F = [ 0 ; f(X(2:N+1),t) ; 0 ];
    
    lado_derecho = U + deltat*F;
    
    % calculamos la nueva U 
    U = matriz \ lado_derecho;
    
    % graficamos la solucion a tiempo t
    figure(1);  
    plot(X,U(1:N+1),'*-');
    title(sprintf('t = %5.3f',t));
    axis(ejes);  pause(0.01);    
end