%第一步：寫一個m文件 brownfgh.m 計算目標函數和目標函數的梯度以及稀疏海色矩陣。
function [f,g,H] = brownfgH(x)
%計算目標函數
  n=length(x); y=zeros(n,1);
  i=1:(n-1);
  y(i)=(x(i).^2).^(x(i+1).^2+1)+(x(i+1).^2).^(x(i).^2+1);
  f=sum(y);
%計算函數梯度矩陣。
  if nargout > 1
     i=1:(n-1); g = zeros(n,1);
     g(i)= 2*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
            2*x(i).*((x(i+1).^2).^(x(i).^2+1)).*log(x(i+1).^2);
     g(i+1)=g(i+1)+2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).* ... 
log(x(i).^2)+2*(x(i).^2+1).*x(i+1).*((x(i+1).^2).^(x(i).^2));
  end
%
% 計算稀疏對稱海色矩陣
if nargout > 2
v=zeros(n,1);
for i=1:(n-1);
v(i)=2*(x(i+1).^2+1).*((x(i).^2).^(x(i+1).^2))+4*(x(i+1).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i).^2).^ ... 
((x(i+1).^2)-1))+ 2*((x(i+1).^2).^(x(i).^2+1)).*(log(x(i+1).^2));
v(i)=v(i)+4*(x(i).^2).*((x(i+1).^2).^(x(i).^2+1)).*((log(x(i+1).^2)).^2);
v(i+1)=v(i+1)+2*(x(i).^2).^(x(i+1).^2+1).*(log(x(i).^2))+.4*(x(i+1).^2).*((x(i).^2).^(x(i+1).^2+1)).* ..
. ((log(x(i).^2)) .^2)+2*(x(i).^2+1).*((x(i+1).^2).^(x(i).^2));  
v(i+1)=v(i+1)+4*(x(i).^2+1).*(x(i+1).^2).*(x(i).^2).* ((x(i+1).^2).^(x(i).^2-1));
v0=v;
v=zeros(n-1,1);
v(i)=4*x(i+1).*x(i).*((x(i).^2).^(x(i+1).^2))+4*x(i+1).*(x(i+1).^2+1).*x(i).*((x(i).^2).^ ... 
(x(i+1).^2)).*log(x(i).^2);
 v(i)=v(i)+4*x(i+1).*x(i).*((x(i+1).^2).^(x(i).^2)).*log(x(i+1).^2);
v(i)=v(i)+4*x(i).*((x(i+1).^2).^(x(i).^2)).*x(i+1);
v1=v;
i=[(1:n)';(1:(n-1))'];
    j=[(1:n)';(2:n)'];
s=[v0;2*v1];
H=sparse(i,j,s,n,n);
H=(H+H')/2;
end
end
%第二步：使用命令fmincon 解決此問題。
fun = @brownfgh;
load browneq % 得到A和B的等式
n = 1000;
xstart = -ones(n,1); xstart(2:2:n) = 1;
options = optimset('GradObj','on','Hessian','on', ... 
'PrecondBandWidth', inf);
[x,fval,exitflag,output] = ...
fmincon(fun,xstart,[],[],Aeq,beq,[],[],[],options);