function aufg22
% exercise 22
% numerical approximation of a reaction diffusion equation 
% Raphael Kruse, 16. Jul 2009

% parameters
t_0     = 0; % time boundary
t_end   = 1;
dt      = [0.01, 0.005, 0.001, 0.0001]; 

x_0     = 0; % spatially boundary
x_end   = 1;
dx      = [0.1 0.1 0.05 0.02];

gamma_0 = 1; % boundary values
gamma_1 = 1.25;

c = 2;

% time step size for plotting
dtp = 0.05;

% initial value
u_0     = @(x) (gamma_1 - gamma_0)*x + gamma_0;

for i = 1:length(dx)
    % generating the line system
    
    N = (x_end - x_0)/dx(i) - 1;
    % matrix A
    A = (N+1)^2*(diag(ones(1,N-1),1) -2*diag(ones(1,N),0)...
        +diag(ones(1,N-1),-1));
    % vector g
    g = zeros(N,1);
    g(1)   = gamma_0;
    g(end) = gamma_1;
    
    % time and space grid
    t_grid{i} = t_0 : dt(i) : t_end;
    x_grid{i} = x_0 : dx(i) : x_end;    
    
    % allocating memory
    u_expl{i} = zeros(N,length(t_grid{i}));
    u_impl{i} = zeros(N,length(t_grid{i}));
    
    % generating initial value
    u_expl{i}(:,1) = u_0(x_grid{i}(2:N+1));
    u_impl{i}(:,1) = u_0(x_grid{i}(2:N+1));
    
    % applying explicit Euler scheme
    for j = 2:length(t_grid{i})
        u_expl{i}(:,j) = u_expl{i}(:,j-1) + dt(i)*( A * u_expl{i}(:,j-1)...
            -c*u_expl{i}(:,j-1) + g*(N+1)^2);
    end
    
    % applying implicit Euler scheme
    [L,U] = lu((1+dt(i)*c)*eye(N) - dt(i)*A); % LU decomposition
    for j = 2:length(t_grid{i})
        u_impl{i}(:,j) = U\(L\(u_impl{i}(:,j-1) + dt(i)*g*(N+1)^2));
    end
    
end

% plotting
for i = 1:length(dx)
    figure(i)
    k = dtp/dt(i);
    [X,Y] = meshgrid(x_grid{i}(2:end-1),t_grid{i}(1:k:end));
    subplot(1,2,1)
    surf(X,Y,u_expl{i}(:,1:k:end)')
    title(['expl. Euler method, \Delta t = ' num2str(dt(i)) ', \Delta x = '...
        num2str(dx(i))])
    xlabel('x')
    ylabel('t')
    zlabel('u(t,x)')
    subplot(1,2,2)
    surf(X,Y,u_impl{i}(:,1:k:end)')
    title(['impl. Euler method, \Delta t = ' num2str(dt(i)) ', \Delta x = '...
        num2str(dx(i))])
    xlabel('x')
    ylabel('t')
    zlabel('u(t,x)')
end

end
     
    
    
    
    