function [x] = rich(a, b, c, u)
    n = size(b, 1)
    e = zeros(n, 1)
    f = zeros(n, 1)
    x = zeros(n, 1)

    tol = 1e-12

    if b(1) < tol then
        msg = sprintf("b(%d) = %f", i , b);             
        error(msg)
    end
    e(1) = -c(1)/b(1)
    f(1) =  u(1)/b(1)

    for i = 2:n
        DD = a(i) * e(i-1) + b(i)
        if (abs(DD) < tol) then
            msg = sprintf("DD(%d) = %f", i , DD);             
            error(msg)
        end
        e(i) = - c(i) / DD; // e(n) pas utilise
        f(i) = (u(i) - a(i) * f(i-1)) / DD
    end

    x(n) = f(n)

    for j = n-1:-1:1
        x(j) = e(j) * x(j+1) + f(j)
    end

endfunction