function [x,v] = init_Het_Het(mapfile, X, P,ap, eps1,eps2,n)

global hetmds cds
% 
[n1,n2]=size(X);
% X=(x1,..., xn), xi,i=1,...,n are  mesh points 
hetmds.nphase=size(X(1),1);
% P=(p1,...,pn), pi are the parameter vectors corresonding to xi,i=,1,...,n
for j=1:n1
    p(i)=P(:,i);
end
for i=1:n1
p=p(i);
x(i) = [X;p(ap)];
end
v0 = [];
for i=1:n1
hetmds.P0(i) =p(i);
end
for i=1:n1
hetmds(i).ActiveParams = ap(i);
end
func_handles = feval(mapfile);
symord = 0; 
symordp = 0;

if     ~isempty(func_handles{9}),   symord = 5; 
elseif ~isempty(func_handles{8}),   symord = 4; 
elseif ~isempty(func_handles{7}),   symord = 3; 
elseif ~isempty(func_handles{5}),   symord = 2; 
elseif ~isempty(func_handles{3}),   symord = 1; 
end
if     ~isempty(func_handles{6}),   symordp = 2; 
elseif ~isempty(func_handles{4}),   symordp = 1; 
end
if isempty(cds)
    cds.options = contset();
end
cds.options = contset(cds.options, 'SymDerivative', symord);
cds.options = contset(cds.options, 'SymDerivativeP', symordp);

hetmds.mapfile = mapfile;
func_handles = feval(hetmds.mapfile);
hetmds.func = func_handles{2};
hetmds.Jacobian  = func_handles{3};
hetmds.JacobianP = func_handles{4};
hetmds.Hessians  = func_handles{5};
hetmds.HessiansP = func_handles{6};
hetmds.Der3      = func_handles{7};
hetmds.Der4      = func_handles{8};
hetmds.Der5      = func_handles{9};
hetmds.Niterations=n;
siz = size(func_handles,2);
if siz > 9
    j=1;
    for i=10:siz
        hetmds.user{j}= func_handles{i};
        j=j+1;
    end
end
%hetmds.v = [];

func_handles = feval(mapfile);
symord = 0; 
symordp = 0;

if     ~isempty(func_handles{9}),   symord = 5; 
elseif ~isempty(func_handles{8}),   symord = 4; 
elseif ~isempty(func_handles{7}),   symord = 3; 
elseif ~isempty(func_handles{5}),   symord = 2; 
elseif ~isempty(func_handles{3}),   symord = 1; 
end
if     ~isempty(func_handles{6}),   symordp = 2; 
elseif ~isempty(func_handles{4}),   symordp = 1; 
end
if isempty(cds)
    cds.options = contset();
end
cds.options = contset(cds.options, 'SymDerivative', symord);
cds.options = contset(cds.options, 'SymDerivativeP', symordp);


homds.mapfile = mapfile;
homds.func = func_handles{2};
homds.Jacobian  = func_handles{3};
homds.JacobianP = func_handles{4};
homds.Hessians  = func_handles{5};
homds.HessiansP = func_handles{6};
homds.Der3 = func_handles{7};
homds.Der4 = func_handles{8};
homds.Der5 = func_handles{9};

hetmds.Niteration=n;
hetmds.eps0=aps0;
hetmds.eps1=eps1;

A1 =  contjac(x(1));
AN =  contjac(x(n1));
D1 = eig(A1);
DN=eig(AN);
% nneg = dimension of stable subspace
hetmds.nunst=0;
for i=1:hetmds.nphase
    if (abs(D1(i))>1)
        hetmds.nunst=hetmds.nunst+1;
    end
end
hetmds.nst=0;
for i=1:hetmds.nphase
    if (abs(DN(i))<1)
        hetmds.nst=hetmds.nst+1;
    end
end

[Q,R]=schur


% COMPOSE X0
% % ----------

% 1. cycle 
x1 = x(1:hetmds.ncoords);

[x1,v]=Hom_new_mesh(x1,v,ntst,ncol);
ups = reshape(x1,hetmds.nphase,hetmds.tps);
hetmds.upold = ups;

% 2. equilibrium coordinates
x1 = [x1; hetmds.x0];
hetmds.PeriodIdx = length(x1);
% 3. (two) free parameters
x1 = [x1; hetmds.P0(ap)];
% 4. extra free parameters
extravec = [hetmds.T; hetmds.eps0; hetmds.eps1];
x1 = [x1; extravec(find(hetmds.extravec))];
% 5. YS and YU, initialized to 0
for i=1:hetmds.nneg
    x1 = [x1; zeros(hetmds.npos,1)];
end
for i=1:hetmds.npos
    x1 = [x1; zeros(hetmds.nneg,1)];
end

x = x1;
v = [];


% ASSIGN SOME VALUES TO HOMOCLINIC FIELDS
% ---------------------------------------

hetmds.YS = zeros(hetmds.npos,hetmds.nneg);
hetmds.YU = zeros(hetmds.nneg,hetmds.npos);

% Third parameter = unstable_flag, 
% 1 if we want the unstable space, 0 if we want the stable one
[QS, eigvlS, dimS] = computeBase(A,0,hetmds.nneg);
[QU, eigvlU, dimU] = computeBase(A,1,hetmds.npos);

hetmds.oldStableQ = QS;
hetmds.oldUnstableQ = QU;
hetmds.ups = [];
hetmds.upold = [];
hetmds.upoldp = [];

%-----------------------------------------------------------------
function init_hetmds(mapfile,x,p,ap,ntst, ncol,extravec,T,eps0,eps1,s)
global hetmds
oldhetmds = hetmds;
hetmds = [];

hetmds.mapfile = mapfile;
func_handles = feval(hetmds.mapfile);
hetmds.func = func_handles{2};
hetmds.Jacobian  = func_handles{3};
hetmds.JacobianP = func_handles{4};
hetmds.Hessians  = func_handles{5};
hetmds.HessiansP = func_handles{6};
hetmds.Der3=[];
siz = size(func_handles,2);
if siz > 9
    j=1;
    for k=10:siz
        hetmds.user{j}= func_handles{k};
        j=j+1;
    end
else hetmds.user=[];end
hetmds.nphase = oldhetmds.nphase;
% hetmds.x0 = oldhetmds.x0;
hetmds.x0 = x(oldhetmds.ncoords+1:oldhetmds.ncoords+oldhetmds.nphase);
hetmds.ActiveParams = ap;
hetmds.P0 = p;
hetmds.extravec = extravec;
Hom_set_ntst_ncol(ntst,ncol,s.data.timemesh);
hetmds.T = T;
hetmds.eps0 = eps0;
hetmds.eps1 = eps1;
hetmds.cols_p1 = 1:(hetmds.ncol+1);
hetmds.cols_p1_coords = 1:(hetmds.ncol+1)*hetmds.nphase;
hetmds.ncol_coord = hetmds.ncol*hetmds.nphase;
hetmds.col_coords = 1:hetmds.ncol*hetmds.nphase;
hetmds.pars = hetmds.ncoords+(1:3);
hetmds.phases = 1:hetmds.nphase;
hetmds.ntstcol = hetmds.ntst*hetmds.ncol;
hetmds.wp = kron(hetmds.wpvec',eye(hetmds.nphase));
hetmds.pwwt = kron(hetmds.wt',eye(hetmds.nphase));
hetmds.pwi = hetmds.wi(ones(1,hetmds.nphase),:);

hetmds.bialt_M1 = [];
hetmds.bialt_M2 = [];
hetmds.bialt_M3 = [];
hetmds.bialt_M4 = [];
hetmds.multipliers = nan;
hetmds.monodromy = [];
hetmds.multi_r1 = [];
hetmds.multi_r2 = [];
hetmds.ups = [];
hetmds.vps = [];
hetmds.tsts = 1:hetmds.ntst;
hetmds.cols = 1:hetmds.ncol;

hetmds.HTPstep = 0;