function j=cjacp(mapfile,x,p,ap)
global cds
nphase=size(ap,1);
if cds.options.SymDerivativeP >= 1
  j = feval(odefile, 0, x, 'jacobianp', p{:});
  j = j(:,ap);
else
  for i=ap
    p1 = p; p1{i} = p1{i}-cds.options.Increment;
    p2 = p; p2{i} = p2{i}+cds.options.Increment;
    j(:,i) = feval(mapfile, 0, x, '', p2{:})-feval(odefile, 0, x, '', p1{:});
  end
  j = j/(2*cds.options.Increment);
  j=j(:,ap);
end