% Programita para resolver la ecuacion de Poisson
% con condiciones de Dirichlet.
% Metodo: Diferencias finitas
%
%  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 = 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); 
        H1 = 0;
        H2 = 0;
        uE = @(t)(0);
        f = @(x,t)(zeros(size(x)));
        u0 = @(x)(sin(pi*x/2));
    case 2
        L = 1;
        T = 3;
        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 = 4;
        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 = 20;
h = L/N;
lambda = 1/4; % la mitad de lo necesario para tener estabilidad
deltat = h^2/(lambda*k); 


%% 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);

for t=deltat:deltat:T
    % calculamos U en los nodos "interiores"
    U(2:N+1) = deltat*f(X(2:N+1), t) ...
        + lambda*U(1:N) + (1-2*lambda)*U(2:N+1) + lambda*U(3:N+2);
    % le damos valor a U en el extremo izquierdo
    U(1) = a(t);
    % le damos valor a la U en el nodo ficticio
    U(N+2) = (H2*uE(t)-H1*U(N+1))*2*h/k + U(N);
    
    % 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)
        pause
    else
        plot(X,U(1:N+1),'*-');
    end
    title(sprintf('t = %5.3f',t))
    %axis(ejes)
    pause(0.02);
end