function y=f(z)
    y(1)=X1f*z(1)+(1-X1f)*z(2)-1
    y(2)=P*z(1)-10^(6.87776-1171.53/(224.366+z(3)))
    y(3)=P*z(2)-10^(6.923746-1355.126/(209.517+z(3)))
endfunction
function bolha()
    A=input("Digite, sob a forma de matriz linha, os valores de [X1f,P]:")
    X1f=A(1,1)
    P=A(1,2)
    To=input("Digite uma estimativa para a temperatura de ponto de bolha:")
    k1o=10^(6.87776-1171.53/(224.366+To))/P
    k2o=10^(6.923746-1355.126/(209.517+To))/P
    zo=[k1o;k2o;To]
    Tbolha=fsolve(zo,f)
    disp("Temperatura de ponto de bolha:",Tbolha(3,1))
    [Tbolha,X1f]=return(Tbolha,X1f)
endfunction
function y1=f1(z1)
    y1(1)=X1f*F-z1(1)*z1(3)-z1(2)*z1(4)
    y1(2)=(1-X1f)*F-z1(1)*z1(5)-z1(2)*z1(6)
    y1(3)=z1(4)+z1(6)-1
    y1(4)=z1(3)+z1(5)-1
    y1(5)=k11*z1(4)-z1(3)
    y1(6)=k21*z1(6)-z1(5)    
endfunction
function balmass()
    A=input("Digite, sob a forma de matriz linha, os valores de [T1,X1f,F,P]:")
    T1=A(1,1)
    X1f=A(1,2)
    F=A(1,3)
    P=A(1,4)
    k11=10^(6.87776-1171.53/(224.366+T1))/P
    k21=10^(6.923746-1355.126/(209.517+T1))/P
    inicial=[0.5*F;0.5*F;0.5;0.5;0.5;0.5]
    result=fsolve(inicial,f1)
    disp("[V1;L1;Y11;X11;Y21;X21]:",result)
    [result,F,T1]=return(result,F,T1)
endfunction
function balener()
    bolha
    balmass
    h1f=51.72*Tbolha(3,1)
    h2f=66.07*Tbolha(3,1)
    Y11=result(3,1)
    H11=7678+31.83*T1+0.0903*T1^2/2
    V1=result(1,1)
    Y21=result(5,1)
    H21=10444.7+41.836*T1+0.1218*T1^2/2
    X11=result(4,1)
    h11=51.72*T1
    L1=result(2,1)
    X21=result(6,1)
    h21=66.07*T1
    Q1=-F*h1f*X1f-F*h2f*(1-X1f)+V1*Y11*H11+V1*Y21*H21+L1*X11*h11+L1*X21*h21
    disp("Q1=",Q1)
    Q1=return(Q1)
endfunction