function [ u ] = tent1d( x, uab, f )
    n = length( x ) - 2;
    idx = 1./diff(x);

    M = diag( -(idx( 1:end - 1 ) + idx( 2: end )) ) + ...
        diag( idx( 2:end - 1 ), -1 ) +                ...
        diag( idx( 2:end - 1 ),  1 );
     
    b = zeros ( n, 1 );
    
    for i = 1:n
        b(i) = int( f, x(i), x(i + 1), x(i + 2) );
    end

    b(1)   = b(1)   -   idx(1).*uab(1);
    b(end) = b(end) - idx(end).*uab(2);
    
    u = [uab(1); M \ b; uab(2)];
end

% Integrate a function using Romberg integration

function intf = int( f, a, b, c )
    n = 16;
    A = zeros( 4, 4 );
    x1 = linspace( a, b, n + 1 );
    x2 = linspace( b, c, n + 1 );
    h1 = (b - a)/n;
    h2 = (c - b)/n;
    fx1 = f( x1 ).*(x1 - a)/(b - a);
    fx2 = f( x2 ).*(x2 - c)/(b - c);
    A(1, 1) = 0.5*(fx1(1) + 2*sum(fx1(2:end - 1)) + fx1(end))*h1 + ...
              0.5*(fx2(1) + 2*sum(fx2(2:end - 1)) + fx2(end))*h2; 
    
    for i = 2:4
        n = 2*n;
        h1 = h1/2;
        h2 = h2/2;
        x1 = linspace( a, b, n + 1 );
        x2 = linspace( b, c, n + 1 );
        fx1 = f( x1 ).*(x1 - a)/(b - a);
        fx2 = f( x2 ).*(x2 - c)/(b - c);
        A(i, 1) = 0.5*(fx1(1) + 2*sum(fx1(2:end - 1)) + fx1(end))*h1 + ...
              0.5*(fx2(1) + 2*sum(fx2(2:end - 1)) + fx2(end))*h2; 
        
        for j = 2:i
            A(i, j) = (2^(j - 1)*A(i, j - 1) - A(i - 1, j - 1))/(2^(j - 1) - 1);
        end
    end
    
    intf = A(4, 4);
end
