Q = 1.2;
g = 9.81;
b = 1.8;
h0 = 6;
H = 0.075;
c0 = 1;
c1 = H - h0 -(Q .^ 2 / (2 * g * (b * h0).^2))
c2 = 0
c3 = (Q .^2 /(2 * g * b .^ 2)
%f(x) - the function of the polynomil
f=@(x) c0 * x .^ 3 + c1 * x .^ 2 + c3;
%function to find the derivative of the polynomial
function d = derivative(x)
    h = 0.000001;
    d = (f(x + h) - f(x)) / h

newton_raphson = @(x) x - (f(x) / derivative(x))

funcion x = iterate(p, n)
    x = p;
    for i in 1:n
        x = newton_raphson(x)
    end
fprintf('Approximate head is %.15f',x)
end
iterate(80, 100)

%{ another solution:
p = [c0 c1 c2 c3];
r = roots(p);
h = r(1);
fprintf('Approximate head is %.15f',h)
%}


x = 0:0.1:10;
plot(x,f(x))
ylim([-10 10])
xlabel('The head over pump')
ylabel('Bernoulli equation')
grid on
hold on