function f = BVP_Het(x,YS,YU,eps0,p)
global hetds 

% extract ups,  [x0,x1,...,xn]
ups = reshape(x,hetds.nphase,hetds.npoints);
p = num2cell(p);

% -----------
% Fixed points condition
J=hetds.niteration;
x0=ups(:,1);
xx=x0;
for i=1:J
xx=feval(hetds.func,0,xx,p{:}); %f^i(x0)-x0
end 
f(1:hetds.nphase,1)=xx-x0;
for i=2:hetds.npoints-2
     xx=ups(:,i);
   for k=1:J
     xx=feval(hetds.func,0,xx,p{:});
   end    
    f(end+1:end+hetds.nphase,1)=xx-ups(:,i+1);%f^i(xi)-xi+1
end

 x1=ups(:,end);xx=x1;
for i=1:J
  xx=feval(hetds.func,0,xx,p{:}); 
end 
  f(end+1:end+hetds.nphase,1)=xx-x1;%f^i(x1)-x1

% -----------
% Ricatti equations
%size(f),
f(end+1:end+hetds.nu,1) = RicattiEval(x0,p,1,YU);
f(end+1:end+hetds.ns,1) = RicattiEval(x1,p,0,YS);

%size(f),pause
% -----------
% Last and first vectors along stable and unstable eigenspaces
Q0U = hetds.Q0;
Q1S = hetds.Q1;
if hetds.nu
    Q1U = Q0U * [-YU'; eye(size(YU,1))];
    
    for i=1:hetds.nu
        f(end+1:end+i,1) = (ups(:,2) - ups(:,1))' * Q1U(:,end-i+1);
    end
end

if hetds.ns
    Q1S = Q1S * [-YS'; eye(size(YS,1))];
    for i=1:hetds.ns
        f(end+1:end+i,1) = (ups(:,end-1) - ups(:,end))' * Q1S(:,end-i+1);
    end
end

% -----------
% Distances from endpoints to fixed points equal to epsilons
eps1=hetds.eps1;
 f(end+1,1) = norm(ups(:,2) - x0)-eps0;
 f(end+1,1) = norm(ups(:,end-1) - x1)-eps1;
  f,NF=norm(f),pause
 %normf=norm(f),pause