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=1e-4*ones(m,1);  % 1e-4 (e nao zero) para aparecer na figura
cs=contour(x1,x2,L,[-5 0 5 10 15 18 19 19.4]);
xlabel('X1');
ylabel('X2');
clabel(cs);
hold on;
plot(x1(1:m),r2,'b:');
plot(r1,x2(1:n),'b:');
p=min(n,m);
plot(x1(1:p),x2(1:p),'b:');
figure(gcf);
hold off;
% cinvert
figure(2); mesh(X1,X2,L);
