function [ u ] = uniform1d( 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) = 0.5*int( f, x(i), 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 )
    n = 16;
    A = zeros( 4, 4 );
    x = linspace( a, b, n + 1 );
    h = (b - a)/n;
    fx = f( x );
    A(1, 1) = 0.5*(fx(1) + 2*sum(fx(2:end - 1)) + fx(end))*h; 
    
    for i = 2:4
        n = 2*n;
        h = h/2;
        x = linspace( a, b, n + 1 );
        fx = f( x );
        A(i, 1) = 0.5*(fx(1) + 2*sum(fx(2:end - 1)) + fx(end))*h; 
        
        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
