Here I specify x and y as well as the contour levels in the contour3 call. Adjust as desired.
%% Main Script.
%% Clear the MATLAB Workspace.
clear
clc
close all
format compact
format long
% Constants.
mu = 0.0121505856; % Mass parameter.
% Setup the mesh for computing the effective potential some range and
% in steps of n.
X = -2:0.01:2;
Y = X;
[ x, y ] = meshgrid(X, Y);
% Compute the effective potential.
U_bar = Effective_Potential(x, y, mu);
% Plot the contours of the effective potential.
figure(1),
surf(x, y, U_bar, 'EdgeColor', 'none'), ax = gca;
ax.CLim = [1.5 5.0]; grid on, hold on,
contour3(x, y, U_bar, 1.5:0.1:5.0, '-k', 'Linewidth', 1.5), hold off,
title('Earth-Moon Pseudo Potential', 'Fontsize', 15),
xlabel('x'), ylabel('y'), zlabel('z'),
axis([-1.5 1.5 -1.5 1.5 1.49 5.0]), view([26, 40]),
colorbar, colormap('turbo'), shading('interp')
%% Asscoiated functions.
function [ U ] = Effective_Potential(x, y, mu)
% Inputs:
% x: A 2D mesh grid of N by N size for the x-coordinate.
% y: A 2D mesh grid of N by N size for the y-coordinate.
% mu: The mass parameter for a CR3BP system.
% Outputs:
% U_bar: 2D mesh grid of N by N size containing the effective
% potential values for the given x and y coordinates.
% Compute the distance of the s/c from the smaller primary.
r_1 = sqrt((x + mu - 1).^2 + y.^2);
% Compute the distance of the s/c from the bigger primary.
r_2 = sqrt((x + mu).^2 + y.^2);
% Compute the effective potential.
U = -(-(1 - mu)./r_2 - mu./r_1 -(1/2)*(x.^2 + y.^2));
end
