% 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) + H1 u(L,t) = H2 uE(t)
%  u(x,0) = u0(x)
%

%% Parametros del problema
caso = 1
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); 
        H1 = 0;
        H2 = 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));
        H1 = 0;
        H2 = 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);
        H1 = 0;
        H2 = 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);
        H1 = 3;
        H2 = 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 = 10;
h = L/N;
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+2) = (H2*uE(0)-H1*U(N+1))*2*h/k + U(N);

%% 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(N+2, N:N+2) = [-1 2*h*H1/k 1];
matriz(1,1:3) = [1 0 0];


%% Resolución
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(1) = 0;
    F(2:N+1) = f(X(2:N+1),t);
    F(N+2) = 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+1),'*',x,u)
    else
        plot(X,U(1:N+1),'*-');
    end
    title(sprintf('t = %5.3f',t))
    axis(ejes)
    pause(0.01);
end