Untitled

x0 = [0.1 0.01 0.01];  

v_sol = fsolve(@serbatoi, x0);

fprintf('v_sol1 = %.4f m^3\n', v_sol(1));
fprintf('v_sol2 = %.4f m^3\n', v_sol(2));
fprintf('v_sol3 = %.4f m^3\n', v_sol(3));

function eqn = serbatoi(x)
    D1 = 0.3;    % [m]
    D3 = 0.4;    % [m]
    epsilon1 = 0.01*D1;
    epsilon3 = 0.01*D3;
    rho = 1000;  % [kg/m^3]
    viscosita = 1.141e-3;  % [Pa*s]
    za = 80;     % [m]
    zb = 30;     % [m]
    L1 = 300;    % [m]
    L3 = 900;    % [m]
    g = 9.81;    % [m/s^2]

    Q = x(1);
    f1 = x(2);
    f3 = x(3);

    % Calcolo Reynolds
    Re1 = (rho*Q/(pi/4*D1^2)*D1/viscosita);
    Re3 = (rho*Q/(pi/4*D3^2)*D3/viscosita);

    eqn(1) = 1/sqrt(f1) + 4*log10(epsilon1/(D1*3.7) + 1.255/(Re1*sqrt(f1)));
    eqn(2) = 1/sqrt(f3) + 4*log10(epsilon3/(D3*3.7) + 1.255/(Re3*sqrt(f3)));
    eqn(3) = za - zb - (4*L1/D1*f1*(Q/(pi/4*D1^2))^2/(2*g)) - (4*L3/D3*f3*(Q/(pi/4*D3^2))^2/(2*g));

    eqn = eqn(:);
end