ODE Event Function not Properly Configured

I am trying to execute two event functions "reaching L2" and "going beyond L3". But everytime I get this error:

Error using odeevents

The ODE option 'Events' must be set to 'on' | 'off'

Error in ode45 (line 140)

odeevents(odeIsFuncHandle,ode,t0,y0,options,varargin);

Error in Lyapunov_Manifolds_Main (line 195)

[t1u, x1u, te1u, xe1u, ie1u] = ode45('CR3BP', 4*t1span, ..

-----------------------------------------------------------------------------------------------------------------

%% Clear MATLAB Workspace.
clear
clc
close all
format compact
format long
%% Constants.
mu = 0.0121505856;            % Mass parameter of the Earth-Moon system.
n = 1;                        % Mean motion of the Moon's orbit around Earth.
L2x = 1.155682165407869;      % x-coordinate for the L2 Lagrange point.
% Set the time periods of L1 and L2 Lyapunov orbits.
T1 = 2.7953;
T2 = 3.4155;
% Set the initial STM.
phi0 = eye(6,6);
%% Integration setup.
% Input setup.
s1_0 = [               0.8189; 
                            0; 
                            0; 
                            0; 
                      0.17454; 
                            0; 
        reshape(phi0, 36, 1)];
s2_0 = [                1.1809; 
                             0; 
                             0; 
                             0; 
                      -0.15587; 
                             0; 
         reshape(phi0, 36, 1)];
t1span = [0 T1];
t2span = [0 T2];
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
options1 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', 'reaching_L2'); 
options2 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', 'going_beyond_L3');
%% Integrate the L_1 and L_2 Lyapunov orbits.
[~, x1] = ode45('CR3BP_STM', t1span, s1_0, options);
[~, x2] = ode45('CR3BP_STM', t2span, s2_0, options);
%% Start plotting the orbits in a figure, later to overlay with the
%  manifolds in this same figure.
plot3(x1(:,1), x1(:,2), x1(:,3), 'r', 'DisplayName', 'L_1 Lyapunov Orbit',...
    'LineWidth', 3);
hold on;
plot3(x2(:,1), x2(:,2), x2(:,3), 'r', 'LineWidth', 3, ...
    'DisplayName', 'L_2 Lyapunov Orbit');
scatter3(1 - mu, 0, 0, 'b*', 'LineWidth', 4, 'DisplayName', 'Moon');
scatter3(-mu, 0, 0, 'b*', 'LineWidth', 10, 'DisplayName', 'Earth');
grid on;
axis equal;
xlabel('x [NDU]');
ylabel('y [NDU]');
zlabel('z [NDU]');
%% Select 100 points along each of the L1 and L2 Lyapunov orbits 
%  respectively.
points = 100;
x1Points = zeros(points, 6);
x2Points = zeros(points,6);
step1 = floor(size(x1, 1)/points);
step2 = floor(size(x2, 1)/points);
for ii = 1:points
    x1Points(ii,:) = x1((ii - 1)*step1 + 1, 1:6);
    x2Points(ii,:) = x2((ii - 1)*step2 + 1, 1:6);
end
%% Initialize plotting control variables for stable and unstable manifolds
%  of the L1 and L2 Lyapunov orbits.
magenta = 0;
yellow = 0;
red = 0;
cyan = 0;
%% Iterate over the 100 points along the L1 and L2 Lyapunov orbits.
%  First 100 iterations go over in the positive direction of perturbations 
%  for both stable and unstable manifolds for L1 orbit and in the negative
%  direction of perturbations for both stable and unstable manifolds for L2
%  orbit.
%  Second 100 iterations go over in the negative direction of perturbations
%  for both stable and unstable manifolds for L1 orbit and in the positive
%  direction of perturbations for both stable and unstable manifolds for L2
%  orbit.
%  In both the 100 iterations, same 100 points along the original orbits are
%  used.
% L_1.
for ii = 1:2*points
    % Set the correct index for initial points to perturb (especially 
    % when ii greater than 100.
    k = ii;
    if ii > points
        k = ii - points;
    end
end
%% Integration setup.
s1_0 = [x1Points(k,:)'; reshape(phi0, 36, 1)];
s2_0 = [x2Points(k,:)'; reshape(phi0, 36, 1)];
%% Integrate the Lyapunov L1 and L2 orbits.
[t1, x1] = ode45('CR3BP_STM', t1span, s1_0, options);
[t2, x2] = ode45('CR3BP_STM', t2span, s2_0, options);
    
%% Compute the monodromy matrix for L1 and L2 orbits.
monodromy1 = reshape(x1(end, 7:42), 6, 6);
monodromy2 = reshape(x2(end, 7:42), 6, 6);
%% Compute the eigenvectors and digonal matrices of eigenvalues for L1 and
%% L2 orbits.
[eigVectors1, eigVDiag1] = eig(monodromy1);
[eigVectors2, eigVDiag2] = eig(monodromy2);
%% Extract eigenvalues from the diagonal matrices of eigenvalues.
eigValues1 = diag(eigVDiag1);
eigValues2 = diag(eigVDiag2);
%% Extract the indices of the unstable and stable eigenvalues.
[~, uIndex1] = max(eigValues1); 
[~, uIndex2] = max(eigValues2);
[~, sIndex1] = min(eigValues1);
[~, sIndex2] = min(eigValues2);
%% Set the perturbation value, approximately 38 kilometers.
d = 0.0001;
%% Normalize the unstable and stable eigenvectors using the position
%  components of the corresponding eigenvectors.
normUEigVec1 = eigVectors1(:, uIndex1)/norm(eigVectors1(1:3, uIndex1));
normUEigVec2 = eigVectors2(:, uIndex2)/norm(eigVectors2(1:3, uIndex2));
normSEigVec1 = eigVectors1(:, sIndex1)/norm(eigVectors1(1:3, sIndex1));
normSEigVec2 = eigVectors1(:, sIndex1)/norm(eigVectors1(1:3, sIndex2));
%% Introduce perturbations into the initial conditions using the unstable
%  and stable eigenvectors and perturbation value.
if ii < points + 1
    x1InitUnstable = x1Points(ii, :)' + d*normUEigVec1;
    x2InitUnstable = x2Points(ii, :)' - d*normUEigVec2;
    x1InitStable = x1Points(ii, :)' + d*normSEigVec1;
    x2InitStable = x2Points(ii, :)' - d*normSEigVec2;
else
    k = ii - points;
    x1InitUnstable = x1Points(k, :)' - d*normUEigVec1;
    x2InitUnstable = x2Points(k, :)' + d*normUEigVec2;
    x1InitStable = x1Points(k, :)' - d*normSEigVec1;
    x2InitStable = x2Points(k, :)' + d*normSEigVec2;
end    
%% Integrate the perturbed initial conditions forward to get the unstable 
%  manifolds.
[t1u, x1u, te1u, xe1u, ie1u] = ode45('CR3BP', 4*t1span, ...
    x1InitUnstable, options1);
[t2u, x2u, te2u, xe2u, ie2u] = ode45('CR3BP', 4*t1span, ...
    x2InitUnstable, options2);
%% For unstable L1 manifolds, plot only those manifolds that do not go to
%  L2 point, but traverse towards the Earth and the interior region, this
%  is for Condition - 1 from the problem statement.
if size(ie1u, 1) == 0
    p1u = plot3(x1u(:,1), x1u(:,2), x1u(:,3), 'm', 'DisplayName', ...
            'L_1 Lyapunov Unstable'); 
        
    % Conditions for suppressing legend entries in the manifolds plot for 
    % every manifold except the first one.
    if magenta == 1
        set(get(get(p1u, 'Annotation'), 'LegendInformation'),...
            'IconDisplayStyle', 'off');
    else
        magenta = 1;
    end
end
% For Condition - 2 of the problem statement, the following logic is used to
% enforce plotting only those unstable manifolds that go beyond Moon, L3
% and cross the x-axis (y=0).
indicesUnstable1 = find(ie2u == 1);
indicesUnstable2 = find(ie2u == 2);
indicesUnstable3 = find(ie2u == 3);
startUnstable = 1;
if size(indicesUnstable2, 1) > 0
    startUnstable = indicesUnstable2(1);
end
indicesUnstable3 = find(ie2u(startUnstable:end) == 3);
if size(indicesUnstable2, 1) > 0
    indexUnstable = find(abs(t2u - te2u(indicesUnstable2(1))) < 0.005);
    if size(indicesUnstable3, 1)>0
        indexUnstable = find(abs(t2u - te2u(startUnstable-1+...
            indicesUnstable3(1))) < 0.008);
    end
    p2u = plot3(x2u(1 : indexUnstable, 1), x2u(1 : indexUnstable, 2), ...
        x2u(1 : indexUnstable,3), 'r', 'DisplayName', 'L_2 Lyapunov Unstable');
    % Conditions for suppressing legend entries in the manifolds plot for 
    % every manifold except the first one.
    if red == 1
        set(get(get(p2u, 'Annotation'), 'LegendInformation'),...
            'IconDisplayStyle', 'off');
    else
        red = 1;
    end
end
% Integrate the perturbed initial conditions backwards for stable
% manifolds.
[t1s, x1s, te1s, xe1s, ie1s] = ode45('CR3BP', -4*t1span, ...
    x1InitStable, options1);
[t2s, x2s, te2s, xe2s, ie2s] = ode45('CR3BP', -4*t2span, ...
    x2InitStable, options2);
% For stable L1 manifolds, plot only those manifolds that do not go to
% L2 point, but traverse towards the Earth and the interior region, this
% is for Condition - 1 from the problem statement.
if size(ie1s, 1) == 0
    p1s = plot3(x1s(:,1), x1s(:,2) ,x1s(:,3), 'y', 'DisplayName', ...
        'L_1 Lyapunov Stable');
    % Conditions for suppressing legend entries in the manifolds plot for 
    % every manifold except the first one.
    if yellow == 1
        set(get(get(p1s, 'Annotation'), 'LegendInformation'),...
            'IconDisplayStyle', 'off');
    else
        yellow = 1;
    end
end
% For Condition - 2 of the problem statement, the following logic is used to
% enforce plotting only those stable manifolds that go beyond Moon, L3
% and cross the x-axis (y=0).
indicesStable1 = find(ie2s == 1);
indicesStable2 = find(ie2s == 2);
start = 1;
if size(indicesStable2, 1) > 0
    start = indicesStable2(1);
end
indicesStable3 = find(ie2s(start : end) == 3);
if size(indicesStable2, 1) > 0
    indexStable = find(abs(t2s - te2s(indicesStable2(1))) < 0.005);
    if size(indicesStable3, 1)>0
        indexStable=find(abs(t2s - te2s(start - 1 + indicesStable3(1))) < 0.005);
    end
    p2s = plot3(x2s(1 : indexStable, 1),x2s(1 : indexStable, 2),...
            x2s(1 : indexStable,3), 'c', 'DisplayName', 'L_2 Lyapunov Stable');
        
    % Conditions for suppressing legend entries in the manifolds plot for 
    % every manifold except the first one.
    if cyan == 1
        set(get(get(p2s, 'Annotation'), 'LegendInformation'),...
            'IconDisplayStyle', 'off');
    else
        cyan = 1;
    end
end
legend('Location', 'northeast', 'FontWeight', 'bold');
fig.Color = 'white';

Answer 1

[Torsten]

Always work with function handles instead of strings:

options1 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', @reaching_L2); 
options2 = odeset('RelTol', 1e-8, 'AbsTol', 1e-8, 'Events', @going_beyond_L3);

[~, x1] = ode45(@CR3BP_STM, t1span, s1_0, options);
[~, x2] = ode45(@CR3BP_STM, t2span, s2_0, options);

[t1, x1] = ode45(@CR3BP_STM, t1span, s1_0, options);
[t2, x2] = ode45(@CR3BP_STM, t2span, s2_0, options);

[t1u, x1u, te1u, xe1u, ie1u] = ode45(@CR3BP, 4*t1span, ...
    x1InitUnstable, options1);
[t2u, x2u, te2u, xe2u, ie2u] = ode45(@CR3BP, 4*t1span, ...
    x2InitUnstable, options2);

[t1s, x1s, te1s, xe1s, ie1s] = ode45(@CR3BP, -4*t1span, ...
    x1InitStable, options1);
[t2s, x2s, te2s, xe2s, ie2s] = ode45(@CR3BP, -4*t2span, ...
    x2InitStable, options2);

If this does not solve your problem, you will have to include the functions to run your code.

[Mark Munson]

Thank you!