Single shooter Matlab code for finding Halo orbits (*not* validated!)

halofunc.m

function state = halofunc(X)
% Parse state vector into variables
x = X(1);
y = X(2);
z = X(3);
vx = X(4);
vy = X(5);
vz = X(6);
Phi = reshape(X(7:end), 6,6);

% Simplication
mu = Constants.mu;
r1 = sqrt((x+mu)^2 + y^2 + z^2);
r2 = sqrt((x+mu-1)^2 + y^2 + z^2);

ax = x + 2*vy - (1-mu)*(x+mu)/r1^3 - mu*(x+mu-1)/r2^3;
ay = y - 2*vx - (1-mu)*y/r1^3 - mu*y/r2^3;
az = -(1-mu)*z/r1^3 - mu*z/r2^3;
% Compute STM
A = A_matrix(X(1:6)', mu);
Phid = A*Phi;

state = [vx; vy; vz; ax; ay; az; reshape(Phid, 6^2, 1)];

hw8.m

clc
clear
close all
format long

X0 = [1.14 0 -0.16 0 -0.223 0;];
% Found?
X0 = [1.140000000000000 0 2.791036381471859 0 -1.282994906583491 0]
it = 0;
tol = 1e-8;
while 1
    fprintf('iteration %d\n', it)
    tspan = [0 50];
    Phi0 = eye(6);
    ic = [X0, reshape(Phi0, 1, 6^2)];
    options = odeset('Events', @periodevent, 'RelTol', 1e-13, 'AbsTol', 1e-13);
    [T,X] = ode45(@(t,X)halofunc(X), tspan, ic, options);
    atXing = X(end, 1:6)
    if atXing(4) < tol && atXing(6) < tol
        disp(atXing)
        disp(X0)
        disp(T(end))
        plotfunc(X,'Rotating', 'HW8')
        break
    else
        % Let's fix this orbit.
        Phi = reshape(X(end, 7:end), 6, 6);
        dXdt = halofunc(X(end, :)');
        deltaXZDot = [-atXing(4); -atXing(6)];
        el1 = [Phi(4, 3) Phi(4, 5); Phi(6, 3) Phi(6, 5)];
        el2 = (1/dXdt(2)) * [dXdt(4);dXdt(6)]*[Phi(2, 3) Phi(2, 5)];
        rhsMat = el1 - el2;
        deltaState = deltaXZDot' * inv(rhsMat);
        %deltaState = inv(rhsMat) * deltaXZDot;
        X0(3) = X0(3) + deltaState(1);
        X0(5) = X0(5) + deltaState(2);
    end
    it = it + 1;
    if it > 10
        disp('INCOMPLETE')
        plotfunc(X,'Rotating', 'HW8')
        break
    end
end

periodevent.m

function [value, isterminal, direction] = periodevent(t, state)
%value = state(2);
value = t - 1.51074761956563;
isterminal = 1;
direction = -1;
end