function XX = xtaylor(x,h,mue,ord)
%
%  The function taylor computes the taylor coefficients of the restricted
%  three body problem. In order to evitate under- or overflow the j-th
%  coefficient is multiplied with h^j.
%
%  input:
% 	x(N,4)		(intval) initial value of solution, N points
% 	h(N)		initial guess of stepsize
% 	mue		(intval) parameter for mass ratio of the primaries
% 	ord		order of Taylor expansion
%
%  output:
% 	XX(N,4,ord+1)	(intval) Taylor coefficients of solution
%
mu1=1-mue;
mu2=mue;
%
% order zero
%
N=size(x,1);
XX(:,:,1)=x;
%
% order one
%
A1(:,1)=x(:,1)+mu2;
A2(:,1)=x(:,1)-mu1;
B1(:,1)=A1(:,1).^2;
B2(:,1)=A2(:,1).^2;
C(:,1)=x(:,2).^2;
R1(:,1)=B1(:,1)+C(:,1);
R2(:,1)=B2(:,1)+C(:,1);
r13(:,1)=mu1*R1(:,1).^(-1.5);
r23(:,1)=mu2*R2(:,1).^(-1.5);
Uj=A1(:,1).*r13(:,1)+A2(:,1).*r23(:,1);
Vj=x(:,2).*(r13(:,1)+r23(:,1));
ff(:,:,1)=[x(:,3:4), x(:,1)+2*x(:,4)-Uj, x(:,2)-2*x(:,3)-Vj];
XX(:,:,2)=(h*ones(1,4)).*ff(:,:,1);
%
% order loop, order = j
%
for j=2:ord,
  A1(:,j)=XX(:,1,j);
  A2(:,j)=XX(:,1,j);
  B1(:,j)=sum(A1(:,1:j).*A1(:,j:-1:1),2);
  B2(:,j)=sum(A2(:,1:j).*A2(:,j:-1:1),2);
  C(:,j)=sum(XX(:,2,1:j).*XX(:,2,j:-1:1),3);
  R1(:,j)=B1(:,j)+C(:,j);
  R2(:,j)=B2(:,j)+C(:,j);
  r13(:,j)=(R1(:,2:j).*r13(:,j-1:-1:1))*((1-j-[1:j-1]/2)')./(R1(:,1)*(j-1));
  r23(:,j)=(R2(:,2:j).*r23(:,j-1:-1:1))*((1-j-[1:j-1]/2)')./(R2(:,1)*(j-1));
  Uj=sum(A1(:,1:j).*r13(:,j:-1:1) + A2(:,1:j).*r23(:,j:-1:1),2);
  Vj=sum(reshape(XX(:,2,1:j),N,j).*(r13(:,j:-1:1)+r23(:,j:-1:1)),2);
  ff(:,[1:2],j)=XX(:,3:4,j);
  ff(:,3,j)=XX(:,1,j)+2*XX(:,4,j)-Uj;
  ff(:,4,j)=XX(:,2,j)-2*XX(:,3,j)-Vj;
  XX(:,:,j+1)=(h*ones(1,4)).*ff(:,:,j)/j;
end;

return

%
% written   Sept/10/02   D.Stoffer
%

