function Dfw = fprime(x,w,mue)
%
%  The function fprime computes the derivative of the vectorfield.
%
%  input:
% 	x(N,4)		(intval) initial value of solution, N points
% 	w(N,16)		(intval) initial value of derivative, N points
% 	mue		(intval) parameter for mass ratio of the primaries
%
%  output:
% 	Dfw(N,16)	(intval) Df(x)*w
%

N=size(x,1);
mu1=1-mue;
mu2=mue;
R1=(x(:,1)+mu2).^2+x(:,2).^2;
R2=(x(:,1)-mu1).^2+x(:,2).^2;
r1=sqrt(R1); r2=sqrt(R2);
a31=1 - mu1*(x(:,2).^2-2*(x(:,1)+mu2).^2)./r1.^5 - ...
	mu2*(x(:,2).^2-2*(x(:,1)-mu1).^2)./r2.^5;
a32=3*x(:,2).*(mu1*(x(:,1)+mu2)./r1.^5 + mu2*(x(:,1)-mu1)./r2.^5);
a42=1 - mu1*((x(:,1)+mu2).^2-2*x(:,2).^2)./r1.^5 - ...
	mu2*((x(:,1)-mu1).^2-2*x(:,2).^2)./r2.^5;
Dfw(:,[4 8 12 16])=(a32*ones(1,4)).*w(:,[1 5 9 13]) + ...
	(a42*ones(1,4)).*w(:,[2 6 10 14]) - 2*w(:,[3 7 11 15]);
Dfw(:,[1 5 9 13])=w(:,[3 7 11 15]);
Dfw(:,[2 6 10 14])=w(:,[4 8 12 16]);
Dfw(:,[3 7 11 15])=(a31*ones(1,4)).*w(:,[1 5 9 13]) + ...
	(a32*ones(1,4)).*w(:,[2 6 10 14]) + 2*w(:,[4 8 12 16]);

return

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