[MatLab, Metodi Numerici] problema fsolve

[telerino]

Buongiorno a tutti, non so se è questa la sezione giusta o magari Ingegneria con l'etichetta Metodi Numerici, nel caso sposterò. Ho un problema con un sistema di due equazioni non lineari in due incognite, ho esplicitato tutti i termini, scritto nel modo migliore possibile le due equazioni ma niente, matlab mi dice che il problema è incondizionato e, ovviamente, non riesce a trovare una soluzione (l'exitflag mi dà valori sia positivi che negativi, a seconda dei valori di tentativo iniziali). Vi scrivo il codice che ho scritto io e se qualcuno ha qualche idea, suggerimento, dritta avrà tutta la mia riconoscenza
Il problema sorge al secondo fsolve (riga 39), vi allego il programma per intero così tutti i parametri sono calcolati

clear all, clc, close all, format long e
P=[0.048756,0.096854,0.39587,0.8773,1.3636,1.851,2.3294,2.8052,3.2871,3.7559,4.1878];
lp=length(P);
global Tt1 Pt1 Tt2  r2 eps1 eps2 mu11 r1 X12 Pt1tot1 T R  vs1 vs2 P1
Ts1=[305.3];
Ps1=[574.5];
rhos1=[1.51];
Ts2=637.7
Ps2=548.6
rhos2=1.158
T=[353.15];
Mw1=[44];
Mw2=80000;
R=[8.314472];
Tt1=(T/Ts1) %T tilde%
vs1=(R*Ts1)/Ps1; %v specifico%
r1=(Mw1/(rhos1*vs1))
vs2=(R*Ts2)/Ps2;
r2=(Mw2/(rhos2*vs2))
eps1=R*Ts1
eps2=R*Ts2
for i=1:lp
Pt1=(P(i)/Ps1);
Pt1tot(i)=Pt1;
x0 = [0.5];
options=optimset('TolFun',1e-12,'TolX',1e-12)
[rhot1,fval1,exitflag]=fsolve(@puro,x0,options)
rhot1tot(i)=rhot1;
end
rhot1tot
mu1=(-r1.*rhot1tot)./Tt1+((r1.*Pt1)./(Tt1.*rhot1tot))+((r1.*log(1-rhot1tot))./rhot1tot)-(r1.*log(1-rhot1tot))+log(rhot1tot)
k12=0;
X12=(eps1/(R*T))+(eps2/(R*T))-((2*sqrt(eps1*eps2))/(R*T))+((2*k12*sqrt(eps1*eps2))/(R*T));
Pt1tot1=Pt1tot(1)
mu11=mu1(1)
P1=P(1)
x0=[1,1]
options=optimset('TolFun',1e-12,'TolX',1e-12)
[sol,fval,exitflag]=fsolve(@sist_nonlin2,x0,options)

la funzione del primo fsolve è questa:

function F = puro(x)
global Pt1 Tt1 r1
F = (x.^2)+Pt1+Tt1*log(1-x)+Tt1*x-(Tt1*x/r1)
end

mentre per il secondo ho provato due funzioni diverse,
una compatta:

function f=sist_nonlin2(x)
global r1 r2 T R X12 vs1 vs2 eps1 eps2 mu11 P1 Tt1 Pt1tot1
T3=((T*R)/(eps1-x(1)*eps1+x(1)*eps2-x(1)*X12*R*T+(x(1).^2)*X12*R*T));
P3=((P1*vs1-x(1)*P1*vs1+x(1)*P1*vs2)/(eps1-x(1)*eps1+x(1)*eps2-x(1)*X12*R*T+(x(1).^2)*X12*R*T));
r3=((x(1)*r1*r2-r1*r2)/(x(1)*r2-x(1)*r1-r2)+((-x(1)*(r2.^2))/(x(1)*r2-x(1)*r1-r2))+r2);
v3=1/x(2);
f=[(x(2)^2)+P3+T3*log(1-x(2))+T3*x(2)-((T3*x(2))/r3),
log(1-x(1))+x(1)-(r1/r2)*x(1)+r1*x(2)*(x(1)^2)*X12-((r1*x(2))/Tt1)+((r1*Pt1tot1*v3)/(Tt1))+((r1*log(1-x(2))*v3))-r1*log(1-x(2))+log(x(2))-mu11];
end

e una esplicita:

function f=sist_nonlin3(x)
global r1 r2 T R X12 vs1 vs2 eps1 eps2 mu11 P1 Tt1 Pt1tot1
f=[(log(1-x(1))+x(1)-(r1/r2)*x(1)+x(2)*(x(1).^2)*r1*X12-((x(2)*r1)/Tt1)+((r1*Pt1tot1)/(x(2)*Tt1))+((r1*log(1-x(2)))/x(2))-r1*log(1-x(2))+log(x(2))-mu11),
((x(2).^2)+((P1*vs1-x(1)*P1*vs1+x(1)*P1*vs2)/(eps1-x(1)*eps1+x(1)*eps2-x(1)*X12*R*T+(x(1).^2)*X12*R*T))+((T*R)/(eps1-x(1)*eps1+x(1)*eps2-x(1)*X12*R*T+(x(1).^2)*X12*R*T))*log(1-x(2))+((T*R)/(eps1-x(1)*eps1+x(1)*eps2-x(1)*X12*R*T+(x(1).^2)*X12*R*T))*x(2)-((((T*R)/(eps1-x(1)*eps1+x(1)*eps2-x(1)*X12*R*T+(x(1).^2)*X12*R*T))*x(2))/((x(1)*r1*r2-r1*r2)/(x(1)*r2-x(1)*r1-r2)+((-x(1)*(r2.^2))/(x(1)*r2-x(1)*r1-r2))+r2)))]
end

il sitema comunque è il seguente:
$\{(x(2)^2+P_3+T_3*ln(1-x(2))+T_3*x(2)-(T_3*x(2))/r_3=0),(ln(1-x(1))+x(1)-r_1/r_2*x(1)+x(2)x(1)^2r_1X_12-(x(2)r_1)/(Tt_1)+(r_1Pt_1v_3)/(Tt_1)+r_1ln(1-x(2))v_3-r_1ln(1-x(2))+ln(x(2))-mu_1=0):}$

i termini sono tutti precedentemente calcolati e $P_3$,$T_3$,$r_3$ e $v_3$ che sono funzione di $x(1)$ e $x(2)$ nella prima funzione li ho definiti, nella seconda esplicitati (le loro relazioni le trovate nella funzione sist_nonlin2).
Vi ringrazio se riuscirete a darmi una mano.

[telerino]

Ho risolto da solo, scrivere tutto il papiello e perderci ancora la testa mi ha aiutato a risolvere, potete chiudere!