function result = aufgabe2()
    N = 40;
    time1 = zeros(N);
    err1 = zeros(N);
    
    time2 = zeros(N);
    err2 = zeros(N);
    
    for m = 2:2:N
        f = ones(m*m,1);
        A = AN(m);
        
        tic;
        [L, R] = LRZerlegung(A);
        x = Loesen(L, R, f);
        time1(m) = toc;
        
        err1(m) = norm(A*x-f);
        printf("m = %d,\tn = %d:\tt=%f,\tnorm=%f\n", m, m*m, time1(m), err1(m));
    end
    
    printf("-----------------------------------------\n");
    
    for m = 2:2:N
        f = ones(m*m,1);
        A = AN(m);
        
        tic;
        [L, R] = LRZerlegungBand(A);
        x = LoesenBand(L, R, f);
        time2(m) = toc;
        
        err2(m) = norm(A*x-f);
        printf("m = %d,\tn = %d:\tt=%f,\tnorm=%f\n", m, m*m, time2(m), err2(m));
    end
    
    result = NaN
end

function result = AN(m)
    B = 4*eye(m)-diag(ones(m-1,1),1)-diag(ones(m-1,1),-1);
    result = [];
    for i = 1:m
        row = [];
        for j=1:m
            if j == i
                add = B;
            elseif j == i-1 || j == i+1
                add = (-1)*eye(m);
            else
                add = zeros(m);
            end
            row = [row add];
        end
        result = [result; row];
    end
    result;
end

function [L, R] = LRZerlegung(A)
    [rows, cols] = size(A);
    if rows != cols
        L = NaN;
        R = NaN;
        return
    end
    
    n = cols;
    
    for i = 1:n
        A(i+1:n, i) = A(i+1:n, i)/A(i, i);
        A(i+1:n, i+1:n) = A(i+1:n, i+1:n) - A(i+1:n, i) * A(i, i+1:n);
    end
    R = triu(A);
    L = A-R+eye(n);
end

function x = Loesen(L, R, b)
    [cols, rows] = size(R);
    x = zeros(cols,1);
    y = zeros(cols,1);
    
    y(1) = b(1)/L(1,1);
    for i = 2:cols
        y(i) = (b(i) - L(i, 1:i-1)*y(1:i-1))/L(i,i);
    end
    
    x(cols) = y(cols)/R(cols,cols);
    for i = (cols-1):-1:1
        x(i) = (y(i) - R(i, i+1:cols)*x(i+1:cols))/R(i, i);
    end
end

function [L, R] = LRZerlegungBand(A)
    % Wir wissen, daß A eine Bandmatrix mit Bandbreite 2m ist. (m=Sqrt(n))
    [rows, cols] = size(A);
    if rows != cols
        L = NaN;
        R = NaN;
        return
    end
    
    n = cols;
    m = half_bw = sqrt(n);
    
    for i = 1:n
        A(i+1:min(i+m,n), i) = A(i+1:min(i+m,n), i)/A(i, i);
        A(i+1:min(i+m,n), i+1:min(i+m,n)) = A(i+1:min(i+m,n), i+1:min(i+m,n)) - A(i+1:min(i+m,n), i) * A(i, i+1:min(i+m,n));
    end
    R = triu(A);
    L = A-R+eye(n);
end

function x = LoesenBand(L, R, b)
    [cols, rows] = size(R);
    m = half_bw = sqrt(cols);
    x = zeros(cols,1);
    y = zeros(cols,1);
    
    y(1) = b(1)/L(1,1);
    for i = 2:cols
        y(i) = (b(i) - L(i, max(i-m,1):i-1)*y(max(i-m,1):i-1))/L(i,i);
    end
    
    x(cols) = y(cols)/R(cols,cols);
    for i = (cols-1):-1:1
        x(i) = (y(i) - R(i, i+1:min(i+m,cols))*x(i+1:min(i+m,cols)))/R(i, i);
    end
end