function [XX,WW] = taylor(x,w,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
% 	w(N,16)		(intval) initial value of derivative
% 	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
% 	WW(N,16,ord+1)	(intval) Taylor coefficients of derivative
%
mu1=1-mue;
mu2=mue;
%
% order zero
%
N=size(x,1);
XX(:,:,1)=x;
WW(:,:,1)=w;
%
% 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);
r15(:,1)=r13(:,1)./R1(:,1);
r25(:,1)=r23(:,1)./R2(:,1);
F(:,1)=A1(:,1).*r15(:,1)+A2(:,1).*r25(:,1);
a32(:,1)=3*x(:,2).*F(:,1);
a31(:,1)=1-(C(:,1)-2*B1(:,1)).*r15(:,1)-(C(:,1)-2*B2(:,1)).*r25(:,1);
a42(:,1)=1-(B1(:,1)-2*C(:,1)).*r15(:,1)-(B2(:,1)-2*C(:,1)).*r25(:,1);
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);
WW(:,[1,5,9,13,2,6,10,14],2)=(h*ones(1,8)).*w(:,[3,7,11,15,4,8,12,16]);
WW(:,[3,7,11,15],2)=(h*ones(1,4)).*((a31(:,1)*ones(1,4)).*w(:,[1,5,9,13]) ...
	+ (a32(:,1)*ones(1,4)).*w(:,[2,6,10,14]) + 2*w(:,[4,8,12,16]));
WW(:,[4,8,12,16],2)=(h*ones(1,4)).*((a32(:,1)*ones(1,4)).*w(:,[1,5,9,13]) ...
	+ (a42(:,1)*ones(1,4)).*w(:,[2,6,10,14]) - 2*w(:,[3,7,11,15]));
%
% 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));
  r15(:,j)=(r13(:,j)-sum(R1(:,2:j).*r15(:,j-1:-1:1),2))./R1(:,1);
  r25(:,j)=(r23(:,j)-sum(R2(:,2:j).*r25(:,j-1:-1:1),2))./R2(:,1);
  F(:,j)=sum(A1(:,1:j).*r15(:,j:-1:1) + A2(:,1:j).*r25(:,j:-1:1),2);
  a32(:,j)=3*sum(reshape(XX(:,2,1:j),N,j).*F(:,j:-1:1),2);
  a31(:,j)=-sum((C(:,1:j)-2*B1(:,1:j)).*r15(:,j:-1:1) + ...
		(C(:,1:j)-2*B2(:,1:j)).*r25(:,j:-1:1),2);
  a42(:,j)=-sum((B1(:,1:j)-2*C(:,1:j)).*r15(:,j:-1:1) + ...
		(B2(:,1:j)-2*C(:,1:j)).*r25(:,j:-1:1),2);
  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;
  WW(:,[1,5,9,13,2,6,10,14],j+1)=(h*ones(1,8)).*WW(:,[3,7,11,15,4,8,12,16],j);
  WW(:,[3,7,11,15],j+1)=...
	  (h*ones(1,4)).*((a31(:,1)*ones(1,4)).*squeeze(WW(:,[1,5,9,13],j))...
	+ (a32(:,1)*ones(1,4)).*squeeze(WW(:,[2,6,10,14],j)) ...
	+ 2*WW(:,[4,8,12,16],j));
  WW(:,[4,8,12,16],j+1)=...
	  (h*ones(1,4)).*((a32(:,1)*ones(1,4)).*squeeze(WW(:,[1,5,9,13],j))...
	+ (a42(:,1)*ones(1,4)).*squeeze(WW(:,[2,6,10,14],j)) ... 
	- 2*WW(:,[3,7,11,15],j));
  ss=zeros(N,8);
  if isintval(x), ss=intval(ss); end;
  for k=2:j, 
    ss(:,1:4)=ss(:,1:4) + (h*ones(1,4)).*((a31(:,k)*ones(1,4))...
	  .*squeeze(WW(:,[1,5,9,13],j+1-k))...
	+    (a32(:,k)*ones(1,4)).*squeeze(WW(:,[2,6,10,14],j+1-k)));
    ss(:,5:8)=ss(:,5:8) + (h*ones(1,4)).*((a32(:,k)*ones(1,4))...
	  .*squeeze(WW(:,[1,5,9,13],j+1-k))...
	+    (a42(:,k)*ones(1,4)).*squeeze(WW(:,[2,6,10,14],j+1-k)));
  end;
  WW(:,[3,7,11,15,4,8,12,16],j+1)=WW(:,[3,7,11,15,4,8,12,16],j+1)+ss;
  WW(:,:,j+1)=WW(:,:,j+1)/j;
end;

return

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

