Somethin like that:
clear
close all
% solving 10 nonlinear equations
% x = [p1 p2 p3 x1 x2 x3 l1 l2 l3 l4]
n = 10; % number of variables
RQ = 10; % available water
fun = @(x)lagrange(x,RQ);
x0 = ones(n,1); % initial values
x = fsolve(fun,x0);
% compute objective function
f = (12 - x(1))*x(1) + (20 - 1.5*x(2))*x(2) + (28 - 2.5*x(3))*x(3) ...
-(3*x(1)^1.3 + 5*x(2)^1.2 + 6*x(3)^1.15);
P = x(1:3);
X = x(4:6);
Lambda = x(7:10);
function [F] = lagrange(x,RQ)
F(1) = 12 - 2*x(1) - 1.3*3*x(1)^0.3 - x(7);
F(2) = 20 - 3*x(2) - 1.2*5*x(2)^0.2 - x(8);
F(3) = 28 - 5*x(3) - 1.15*6*x(3)^0.15 - x(9);
F(4) = x(7)*0.9*0.4*x(4)^-0.1 - x(10);
F(5) = x(8)*0.8*0.5*x(5)^-0.2 - x(10);
F(6) = x(9)*0.7*0.6*x(6)^-0.3 - x(10);
F(7) = x(1) - 0.4*x(4)^0.9;
F(8) = x(2) - 0.5*x(5)^0.8;
F(9) = x(3) - 0.6*x(6)^0.7;
F(10) = x(4) + x(5) + x(6) - RQ;
end
