Ps = 0.4;// R$/kg-soluto
Px = 0.01;// R$/kg-solvente
F = 10000;// kg-agua/h
xo = 0.02;// kg-soluto/kg-agua
//s=7e-4;   % kg-solvente/kg-agua
s = 0;// kg-solvente/kg-agua
m = 4;// kg-agua/kg-solvente

a = (Ps*F)*xo+((2*Px)*F)/m-(Px*F)*s;
b = ((Px*F)*xo)/m;
c = Ps*F;
d = (Px*F)/m;

//a=130;
//b=0.5;
//c=4000;
//d=25;

x1 = [0.0001:0.0001:0.035]';
x2 = [0.0001:0.0001:0.025]';
//x1=[0.01:0.001:0.02]'';
//x2=[0.006:0.001:0.015]'';

[X1,X2] = meshgrid(x1',x2');
[n,m] = size(X1);
L = a-b./X1-c*X2-d*X1./X2;
r1 = 0.02*ones(n,1);
r2 = 0.0001*ones(m,1);// 1e-4 (e nao zero) para aparecer na figura
contour(x1,x2,L',[-5,0,5,10,15,18,19,19.4]);
//cs=contour(x1,x2,L,[-5 0 1 2 3 4 5]);

xlabel("X1");
ylabel("X2");

set(gca(),"auto_clear","off");
plot(x1(1:m),r2,"b:");
plot(r1,x2(1:n),"b:");
p = min(n,m);
plot(x1(1:p),x2(1:p),"b:");
set(gca(),"auto_clear","on");

// cinvert
//figure(2); mesh(X1,X2,L);

[X1s,X2s] = meshgrid(0.005:0.0001:0.022,0.005:0.0001:0.014);
Ls = a-b./X1s-c*X2s-d*X1s./X2s;
figure(3);
surf(X1s,X2s,Ls)
xlabel("X1");
ylabel("X2");
