F=21095;
K=[66.7,50,0.015,33.3,100];
z=[0.1741,0.5234,0.1574,0.1038,0.0413];
psi0=0.1;
%METODO FUERZA BRUTA
for V=0:0.5:F
L=F-V;
psi=V/F;
for i=1:5
   x(i)=z(i)/(1+psi*(K(i)-1));
   y(i)=K(i)*x(i);
end
if abs(sum(x)-sum(y))<0.0001 
   sprintf('%s','Fuerza bruta')
   psi,L,V,x,y,sum(x),sum(y)
   break;
end
end
%METODO DE RESOLUCION DE NEWTON-RAPHSON
psi=0.3;
for i=1:5
   x(i)=z(i)/(1+psi*(K(i)-1));
   y(i)=K(i)*x(i);
end
while abs(sum(x)-1)>0.0001 
   for i=1:5
   fx(i)=z(i)*(1-K(i))/(1+psi*(K(i)-1));   
   dfx(i)=-(z(i)*(1-K(i))*(K(i)-1))/(1+psi*(K(i)-1))^2;
  	end
	psi=psi-(sum(fx)/(sum(dfx)));
	for i=1:5
   	x(i)=z(i)/(1+psi*(K(i)-1));
   	y(i)=K(i)*x(i);
	end
end
V=psi*F;
L=F-V;
sprintf('%s','Newton-Raphson')
psi,V,L,x,y,sum(x),sum(y)
%METODO DE MINIMOS CUADRADOS GAUSS_NEWTON
psi=0.2;
psi=fsolve('equi',psi0,[],K,z);
V=psi*F;
L=F-V;
for i=1:5
   x(i)=z(i)/(1+psi*(K(i)-1));
   y(i)=K(i)*x(i);
end
sprintf('%s','Gauss-Newton')
psi,L,V,x,y,sum(x),sum(y)
