J3Dynamics.m

function [sys] = J3Dynamics(t, X, params)

sys = zeros(6,1);

mu = params(1);
J2 = params(2);
R = params(3);
J3 = params(4);

sys(1:3) = X(4:6);

x = X(1);
y = X(2);
z = X(3);

% sys(4) = (-mu*x*(2*(x^2+y^2)^3 + 15*J3*R^3*(x^2+y^2)*z + 6*(x^2+y^2)*z^2 ... 
%     - 20*J3*R^3*z^3 + 6*(x^2+y^2)*z^4 + 2*z^6 + 3*J2*R^2*(x^2+y^2-4*z^2)*(x^2+y^2+z^2))) ...
%     /(2*(x^2+y^2+z^2)^(9/2)); %x acceleration
% 
% sys(5) = (-mu*y*(2*(x^2+y^2)^3 + 15*J3*R^3*(x^2+y^2)*z + 6*(x^2+y^2)*z^2 ... 
%     - 20*J3*R^3*z^3 + 6*(x^2+y^2)*z^4 + 2*z^6 + 3*J2*R^2*(x^2+y^2-4*z^2)*(x^2+y^2+z^2))) ...
%     /(2*(x^2+y^2+z^2)^(9/2)); %y acceleration
% 
% sys(6) = (-mu*(J3*R^3*(-3*(x^2+y^2)^2 + 24*(x^2+y^2)*z^2 - 8*z^4) + ...
%     z*(x^2+y^2+z^2)*(3*J2*R^2*(3*(x^2+y^2) - 2*z^2) + 2*(x^2+y^2+z^2)^2))) ...
%     /(2*(x^2+y^2+z^2)^(9/2)); %z acceleration

sys(4:6) = [(-1/2).*x.*(x.^2+y.^2+z.^2).^(-9/2).*(2.*(x.^2+y.^2).^3+15.*J3.* ...
  R.^3.*(x.^2+y.^2).*z+6.*(x.^2+y.^2).^2.*z.^2+(-20).*J3.*R.^3.* ...
 z.^3+6.*(x.^2+y.^2).*z.^4+2.*z.^6+3.*J2.*R.^2.*(x.^2+y.^2+(-4).* ...
  z.^2).*(x.^2+y.^2+z.^2)).*mu,...
  (-1/2).*y.*(x.^2+y.^2+z.^2).^(-9/2).*( ...
  2.*(x.^2+y.^2).^3+15.*J3.*R.^3.*(x.^2+y.^2).*z+6.*(x.^2+y.^2).^2.* ...
z.^2+(-20).*J3.*R.^3.*z.^3+6.*(x.^2+y.^2).*z.^4+2.*z.^6+3.*J2.* ...
R.^2.*(x.^2+y.^2+(-4).*z.^2).*(x.^2+y.^2+z.^2)).*mu,...
(-1/2).*(x.^2+ ...
  y.^2+z.^2).^(-9/2).*(J3.*R.^3.*((-3).*(x.^2+y.^2).^2+24.*(x.^2+ ...
  y.^2).*z.^2+(-8).*z.^4)+z.*(x.^2+y.^2+z.^2).*(3.*J2.*R.^2.*(3.*( ...
  x.^2+y.^2)+(-2).*z.^2)+2.*(x.^2+y.^2+z.^2).^2)).*mu]';

end