function  [J,grad,hess,y]  = objective(u,d)


k1=d(1);
k2=d(2);
cain=2;
cbin=1.5;
V=500;
R=[1 0;0 1];
Fa=u(1);
Fb=u(2);


x0=[0.0578
    0.0551
    0.7825];
opts=optimset('Display','off','TolX',1e-16,'TolFun',1e-16);
x=fsolve(@equations,x0,opts);
ca=x(1);
cb=x(2);
cc=x(3);

%cost function we want to maximize
J=(Fa+Fb)*cc*(Fa+Fb)*cc/(Fa*cain)-0.5*[Fa Fb]*R*[Fa;Fb];
%solving problem as a minimization problem
J=-J;
y=[ca;cb;cc];

if nargout>=2
    [grad,hess]=gradient;
 else
    grad=[];
    hess=[];
    
end

%model equations at steady sate
    function sys=equations(x)
        
        ca=x(1);
        cb=x(2);
        cc=x(3);
        
        sys(1)=Fa*cain-(Fa+Fb)*ca-k1*ca*cb*V;
        sys(2)=Fb*cbin-(Fa+Fb)*cb-k1*ca*cb*V-2*k2*cb^2*V;
        sys(3)=-(Fa+Fb)*cc+k1*ca*cb*V;
        
    end

%something is wrong here
    function [grad,hess]=gradient
        
%       
%          grad=[ Fas - (0.20972911231316554436787669146725*(0.41945822462633108873575338293449*Fas + 0.41945822462633108873575338293449*Fbs))/Fas - (0.08797260110333681431323083566615*(Fas + Fbs))/Fas + (0.20972911231316554436787669146725*(0.41945822462633108873575338293449*Fas + 0.41945822462633108873575338293449*Fbs)*(Fas + Fbs))/Fas^2,
%                Fbs - (0.20972911231316554436787669146725*(0.41945822462633108873575338293449*Fas + 0.41945822462633108873575338293449*Fbs))/Fas - (0.08797260110333681431323083566615*(Fas + Fbs))/Fas];
        
                syms cas cbs ccs
                syms Fas Fbs
                syms k1s k2s
                cas=cas;
                cbs=cbs;
                ccs=ccs;
                Fas=Fas;
                Fbs=Fbs;
                k1s=k1s;
                k2s=k2s;
                
        
                f1=Fas*cain-(Fas+Fbs)*cas-k1s*cas*cbs*V;
                f2=Fbs*cbin-(Fas+Fbs)*cbs-k1s*cas*cbs*V-2*k2s*cbs^2*V;
                f3=-(Fas+Fbs)*ccs+k1s*cas*cbs*V;
        
                f=[f1;f2;f3];
                
                x=[cas;cbs;ccs];
                u=[Fas;Fbs];
                obj=(Fas+Fbs)*ccs*(Fas+Fbs)*ccs/(Fas*cain)-0.5*[Fas Fbs]*R*[Fas;Fbs];
                obj=-obj;
                Jx=jacobian(obj,x);
                Ju=jacobian(obj,u);
                Fx=jacobian(f,x);
                Fu=jacobian(f,u);
        
                grad=(Ju-Jx*inv(Fx)*Fu);
                hess=jacobian(grad,u);
                  Fas=Fa;
                  Fbs=Fb;
                  ccs=cc;
                  cbs=cb;
                  cas=ca;
                  k1s=k1;
                  k2s=k2;
                grad=subs(grad);
                hess=subs(hess);
%                 pause
%         
        
        
        
    end

grad=subs(grad);
hess=subs(hess);


end

