function soln()

x0 = zeros(3, 1);
lb = zeros(3, 1);
options = optimoptions('fmincon', 'display', 'iter-detailed',...
    'SpecifyObjectiveGradient', true,...
    'SpecifyConstraintGradient', true);

sigmag = 0;
k = 0;

[xopt, fopt] = fmincon(@obj, x0, [], [], [], [], lb, [], @con, options)

sigmax = [0.01, 0.1, 0.05];
sigmag = zeros(2, 1);
[~, ~, ~, dg] = prob(xopt);
nc = 2;
nx = 3;

for i = 1:nc
    for j = 1:nx
        sigmag(i) = sigmag(i) + (dg(i, j)*sigmax(j))^2;
    end
end
sigmag = sqrt(sigmag);
k = norminv(.995);

[xopt, fopt] = fmincon(@obj, xopt, [], [], [], [], lb, [], @con, options)


function [f, df] = obj(x)
   [f, ~, df, ~] = prob(x); 
end

function [g, h, dg, dh] = con(x)
    [~, g, ~, dg] = prob(x);
    g = g + k*sigmag;
    dg = dg';
    h = [];
    dh = [];
end


end