Solving symbolic equation set gives wrong answers when compared to MATLAB code

edit

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

https://cloud.sagemath.com/projects/0...

also at

http://preview.tinyurl.com/SageHWhelp

Thank you in advance for the help! This is driving me crazy!

EDIT: I guess I should add some more explanatory text here about the sage worksheet. The answers in that are all correct up until the very last cell that contains the following:

f1(x1,x2,x3)= (yC1^vC1*yD1^vD1) / (yA1^vA1)*(P/P0)^(vC1+vD1-vA1)-K1
f2(x1,x2,x3)= (yC2^vC2) / (yA2^vA2*yB2^vB2)*(P/P0)^(vC2-vB2-vA2)-K2
f3(x1,x2,x3)= (yC3^vC3) / (yA3^vA3*yB3^vB3)*(P/P0)^(vC3-vB3-vA3)-K3
solutions = solve([f1, f2,f3],x1,x2,x3,solution_dict=True)

EDIT again: Wow that is embarrassing, I put up the wrong MATLAB code for comparison. It also appears that I have lost the MATLAB code to compare with. I guess the question is still there, though: Is sage not capable of solving these equations? I can tell you that I compared each one of the final equations f1 through f3 and they are all correct going into the solver as far as values go. It's just that the solver spits back different values than I expect. I had no less than 4 others check the numbers going in, including my professor, although I will say that none of us are very familiar with sage or python.

EDIT again again: I found the MATLAB code. I really need to reorganize my files.

function HW3_d()
%% Embedded Function

function F = HW3_d(x)


%% Given Conditions
P0=100;     % Ambient pressure (kPa)
T0=298.15;  % Ambient temperature (K)
P=500;      % Products pressure (kPa)
T=1600;     % Products temperature (K)


%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
nO2=25+x(1)-x(2)-2*x(3); 
nN2=141-x(2)-x(3); 
nH2O=9; 
nCO2=8-2*x(1); 
% constants are initial concentrations


%% Dissociation Reactions 

% Reaction #1: 2CO2 -> 2CO + O2
% Reaction #2: N2 + 02 -> 2NO
% Reaction #3: N2 + 202 -> 2NO2

% Number of Moles
nCO = 2*x(1);
nNO = 2*x(2);
nNO2 = 2*x(3);
ntot = nCO2 + nH2O + nN2 + nO2 + nCO + nNO + nNO2;


%% Equilibrium constants at 1600 K from Table A.11
K1=exp(-21.656); K2=exp(-10.547); K3=exp(-20.126);


%% Mole Fractions
yA1=nCO2/ntot; yB1=0; yC1=nCO/ntot; yD1=nO2/ntot;
yA2=nN2/ntot; yB2=nO2/ntot; yC2=nNO/ntot; yD2=0;
yA3=nN2/ntot; yB3=nO2/ntot; yC3=nNO2/ntot; yD3=0;


%% Stoichiometric Coefficients
vA1=2; vB1=0; vC1=2; vD1=1;
vA2=1; vB2=1; vC2=2; vD2=0;
vA3=1; vB3=2; vC3=2; vD3=0;


%% Equilibrium Constants
F=[(yC1^vC1*yD1^vD1)/(yA1^vA1*yB1^vB1)*(P/P0)^(vC1+vD1-vA1-vB1)-K1;
    (yC2^vC2*yD2^vD2)/(yA2^vA2*yB2^vB2)*(P/P0)^(vC2+vD2-vA2-vB2)-K2;
    (yC3^vC3*yD3^vD3)/(yA3^vA3*yB3^vB3)*(P ...