function [af]=makeAgeEns(aD,eD,N)

% Check initial monotonicity
if ~all(diff(aD)>=0)
   error('Initial times are not monotonic')
end


af=zeros(length(aD),N);

i=2;j=1;k=1;
while j<N+1

   % Initialize the time series
   aN=(aD(1)+round(eD(1)*randn(1))*ones(length(aD),1));

   % Keep searching for a good next date
   while i<=length(aD)

      p0=aD(i)+round(eD(i)*randn(1000,1));
      pN=p0(p0>aN(i-1));

      % Check if this value is good
      if ~isempty(pN)
	 aN(i)=pN(1); % grab first random value
	 if i==length(aD)
	    i=2; k=1; % if done, reset counters
	    %rflag=true;
	    break;
	 end
	 i=i+1;
	 continue;
      elseif i==2 % Try a new start 
	 aN=(aD(1)+round(eD(1)*randn(1))*ones(length(aD),1));
      else % Back-track and try again
	 i=max(2,i-2);
         %continue;
      end

      % Break out and start over if can't find a good next value
      if k>1e6 
         disp(['prxy t = ' num2str(i)])
         disp(['ens = ' num2str(j)])
         %aN([i-2:i+1])
         error('Too many tries => something wrong with proxy ages')
      end
      k=k+1;
   end
  
   % Check monotonicity
   if ~all(diff(aN)>=0)
      disp(num2str(j))
      error('Code not creating monotonic proxies')
   end
   
   af(:,j)=aN;
   j=j+1;

end


% Different matrix strategy
%aN=bsxfun(@plus,aD,round(eD.*randn(length(eD),1000)));
%find(all(diff(aN,1,1)>=0))

% Could also try just sorting the columns above and replacing problematic values


