% Programita para resolver la ecuacion de Difusion
% con condiciones de Dirichlet y Robin.
% Metodo: Diferencias finitas
% IMPLICITO
%
%  u_t(x,t) - k u_xx(x,t) = f(x,t),   0 < x < L
%  u(0,t) = a(t)
%  k u'(L,t) + kappa1 u(L,t) = kappa2 uE(t)
%  u(x,0) = u0(x)
%

%% Parametros del problema
caso = 4
switch(caso)
    case 1
        % problema de Neumann con solucion exacta
        % u(x,t) = sin(pi*x/2) * exp(-pi^2/4*k*t)
        L = 1;
        T = 3;
        k = 1;
        a = @(t)(0); 
        kappa1 = 0;
        kappa2 = 0;
        uE = @(t)(0);
        f = @(x,t)(zeros(size(x)));
        u0 = @(x)(sin(pi*x/2));
    case 2
        L = 1;
        T = 6;
        k = 1;
        a = @(t)(sin(2*t));
        kappa1 = 0;
        kappa2 = 1;
        uE = @(t)(0);
        f = @(x,t)(zeros(size(x)));
        u0 = @(x)(zeros(size(x)));
    case 3
        L = 1;
        T = 3;
        k = 1;
        a = @(t)(0);
        kappa1 = 0;
        kappa2 = 1;
        uE = @(t)(1);
        f = @(x,t)(zeros(size(x)));
        u0 = @(x)(zeros(size(x)));
    case 4
        L = 1;
        T = 3;
        k = 1;
        a = @(t)(0);
        kappa1 = 3;
        kappa2 = 3;
        uE = @(t)(0);
        f = @(x,t)(10*exp(-100*(x-.5).^2)*(t<1.5));
        u0 = @(x)(zeros(size(x)));
end

%% Parametros del metodo de resolucion
N = 20;
h = L/(N+1);
deltat = 4*h^2/k; % ocho veces lo necesario para tener estabilidad
                  % en el metodo explicito

lambda = deltat*k/h^2; 

%% Condicion inicial
X = [0:h:L]';
U = u0(X);

figure(1)
plot(X,U,'*-')
ejes = axis;
pause

% le damos valor a la U en el nodo ficticio
U(N+3) = (kappa2*uE(0)-kappa1*U(N+1))*2*h/k + U(N+1);

%% Armado de la matriz

unos = ones(N+3,1);
columnas = [-lambda*unos (1+2*lambda)*unos -lambda*unos];

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


%% Resolución
F=zeros(N+3,1);

for t=deltat:deltat:T
    % Ensamblado del lado derecho
    U(1) = a(t);
    U(N+3) = 2*kappa2*h/k*uE(t);
    
    F(1) = 0;
    F(2:N+2) = f(X(2:N+2),t);
    F(N+3) = 0;
    
    lado_derecho = U + deltat*F;
    
    % calculamos U en los nodos "interiores"
    U = matriz \ lado_derecho;
    
    % solucion exacta (solo para el caso 1)
    if (caso==1)
        x = linspace(0,L,3*N);
        u = exp(-k*pi^2/4*t)*sin(pi*x/2);
        plot(X,U(1:N+2),'*',x,u)
    else
        plot(X,U(1:N+2),'*-');
    end
    title(sprintf('t = %5.3f',t))
    axis(ejes)
    pause(0.01);
end