Water Thermodynamic: pressure calculation

Hi everyone, I'm writing the code of the water's pressure.

I wrote the matrix of the constant:

 A= [ A11  A12  A13  A14  A15  A16  A17;
      A21  A22  A23  A24  A25  A26  A27;
      A31  A32  A33  A34  A35  A36  A37;
      A41  A42  A43  A44  A45  A46  A47;
      A51  A52  A53  A54  A55  A56  A57;
      A61  A62  A63  A64  A65  A66  A67;
      A71  A72  A73  A74  A75  A76  A77;
      A81  A82  A83  A84  A85  A86  A87;
      A91  A92  A93  A94  A95  A96  A97;
      A101 A102 A103 A104 A104 A106 A107; ];

I want to write the same C++ code in Matlab:

 double water::C(int i)
 {
     double tau = Ta/T;
     return (i == 0 ? R*T : R*T*(tau - tauc)*pow(tau - taua[i],i-1));
 }
 
 double water::H(int j)
 {
     double factor, sum, rho_aj;
     rho_aj = (j == 0 ? Roa1 : Roaj);
     sum = 0.0;
     factor = Rho - rho_aj;
     for (int i=6; i>0; i--) {
         sum += (A[i][j] + Rho*(i+1)*A[i+1][j]);
         sum *= factor;
     }
     sum += (A[0][j] + Rho*A[1][j]);
     sum += (exp(-E*Rho)*((1.0 - Rho*E)*A[8][j]
                          + Rho*(2.0 - Rho*E)*A[9][j]));
     sum += A[7][j]*pow(factor,7);
     return Rho*Rho*sum;
 }
 
 double water::Pp()
 {
     double P = Rho*R*T;
     for (int i=0; i<7; i++) {
         P += C(i)*H(i);
     }
     return P;
 }

These are the 3 function: the last one is the pressure.

I wrote the H function like this:

 function [H] = H_fun(rho, j)
    
    global R A11 A21 A31 A41 A51 A61 A71 A81 A91 A101 A12 A22 A32 A42 A52 ...
    A62 A72 A82 A92 A102 A13 A23 A33 A43 A53 A63 A73 A83 A93 A103 A14 ...
    A24 A34 A44 A54 A64 A74 A84 A94 A104 A15 A25 A35 A45 A55 A65 A75 ...
    A85 A95 A105 A16 A26 A36 A46 A56 A66 A76 A86 A96 A106 A17 A27 A37 ...
    A47 A57 A67 A77 A87 A97 A107 E Ta rhoA1 rhoAj F1 F2 F3 F4 F5 F6 ...
    F7 F8 a Tp D1 D2 D3 D4 D5 D6 D7 D8 G1 G2 G3 G4 G5 G6 ...
    u0 s0 M Tc pc rhoC T0 tauC



A = [ A11  A12  A13  A14  A15  A16  A17;
      A21  A22  A23  A24  A25  A26  A27;
      A31  A32  A33  A34  A35  A36  A37;
      A41  A42  A43  A44  A45  A46  A47;
      A51  A52  A53  A54  A55  A56  A57;
      A61  A62  A63  A64  A65  A66  A67;
      A71  A72  A73  A74  A75  A76  A77;
      A81  A82  A83  A84  A85  A86  A87;
      A91  A92  A93  A94  A95  A96  A97;
      A101 A102 A103 A104 A104 A106 A107; ];
  

if j==1
    rho_Aj = rhoA1;
else
    rho_Aj = rhoAj;
end

sum = 0.0;
factor = rho - rho_Aj;

for i = [ 7 6 5 4 3 2 1]
    sum = sum + (A(i,j) + rho*(i+1)*A((i+1),j));
    sum = sum*factor;
end

sum = sum + (A(1,j) + rho*A(2,j));
sum = sum + ( exp(-E*rho)*((1.0 - rho*E)*A(9,j) + rho*(2.0 - rho*E)*A(10,j)) );
sum = sum + A(8,j)*(factor^(7));

H = rho*rho*sum;

end

The C function like this:

function [C]= C_fun(T, j)


    global R A11 A21 A31 A41 A51 A61 A71 A81 A91 A101 A12 A22 A32 A42 A52 ...
    A62 A72 A82 A92 A102 A13 A23 A33 A43 A53 A63 A73 A83 A93 A103 A14 ...
    A24 A34 A44 A54 A64 A74 A84 A94 A104 A15 A25 A35 A45 A55 A65 A75 ...
    A85 A95 A105 A16 A26 A36 A46 A56 A66 A76 A86 A96 A106 A17 A27 A37 ...
    A47 A57 A67 A77 A87 A97 A107 E Ta rhoA1 rhoAj F1 F2 F3 F4 F5 F6 ...
    F7 F8 a Tp D1 D2 D3 D4 D5 D6 D7 D8 G1 G2 G3 G4 G5 G6 ...
    u0 s0 M Tc pc rhoC T0 tauC


tauAj = [ 1.544912 2.5 2.5 2.5 2.5 2.5 2.5 ];

tau = Ta / T;

if j==1
    C = R*T;
else
    C = R*T * (tau - tauC)*((tau - tauAj(j))^(j-2));
end

and the Pressure function:

function [P] = P_fun(rho, T)
    
    global R A11 A21 A31 A41 A51 A61 A71 A81 A91 A101 A12 A22 A32 A42 A52 ...
    A62 A72 A82 A92 A102 A13 A23 A33 A43 A53 A63 A73 A83 A93 A103 A14 ...
    A24 A34 A44 A54 A64 A74 A84 A94 A104 A15 A25 A35 A45 A55 A65 A75 ...
    A85 A95 A105 A16 A26 A36 A46 A56 A66 A76 A86 A96 A106 A17 A27 A37 ...
    A47 A57 A67 A77 A87 A97 A107 E Ta rhoA1 rhoAj F1 F2 F3 F4 F5 F6 ...
    F7 F8 a Tp D1 D2 D3 D4 D5 D6 D7 D8 G1 G2 G3 G4 G5 G6 ...
    u0 s0 M Tc pc rhoC T0 tauC

A = [ A11  A12  A13  A14  A15  A16  A17;
      A21  A22  A23  A24  A25  A26  A27;
      A31  A32  A33  A34  A35  A36  A37;
      A41  A42  A43  A44  A45  A46  A47;
      A51  A52  A53  A54  A55  A56  A57;
      A61  A62  A63  A64  A65  A66  A67;
      A71  A72  A73  A74  A75  A76  A77;
      A81  A82  A83  A84  A85  A86  A87;
      A91  A92  A93  A94  A95  A96  A97;
      A101 A102 A103 A104 A104 A106 A107; ];

    
tauAj = [ 1.544912 2.5 2.5 2.5 2.5 2.5 2.5 ];
      
      
P = rho*R*T;

for j = 1:7
    [C(j)] = C_fun(T,j);
    [H(j)] = H_fun(rho,j);
    P = P + C(j)*H(j);
end
    
end

If I put the right values of rho and T and start the programme, it doesn't return the right pressure (as I see in the water tables)

I think I've made a mistake in the index i,j, when I transform the C++ in code into Matlab code.

Is there someone who can help me? Thank you