Writing sum of two functions for ODE solver

I am trying to write a function to feed to and ode solver. Problem is that the function contains the sum of some variables and another function I computed symbolically and I am not sure about its syntax. Could you please help me understand? Below is the code which runs all the way until it tries to perform said sum. The problem occurs in function SCRTBP. Thank you!

%% Simulation parameters

% System parameters

G = 6.67259e-20;

M_sys = 5.28e11; %[kg] mass of the system

mass_ratio = 0.0093; % mass ratio

mu = 1 / ((1/mass_ratio) + 1);

L_sys = 1.2; % [km] source https://dart.jhuapl.edu/Mission/index.php#

OM_sys = sqrt(G * M_sys / L_sys^3);

% Didymain parameters

m_main = (1 - mu) * M_sys;

C = [-0.016 1.82e-2 5.5e-3 -2.06e-5 5.91e-7];

R_main = 780e-3;

% Didymoon parameters

m_moon = mu * M_sys;

a_moon = 103e-3;

b_moon = 79e-3;

c_moon = 66e-3;

rho_moon = 2300;

% Settings

options = odeset('Reltol',1.e-13,'Abstol',1.e-20);

%% SCR3BP

l = 2;

main_p = [G m_main R_main l];

moon_p = [a_moon b_moon c_moon G rho_moon];

sys_p = [OM_sys L_sys mu];

xx0L1_ok = [0.844746604107861 0 0.0825903835426480 0 0.170006967699603 -0.00962175624416549];

TL1_ok = 2.76472515631288;

sol01s = ode113(@(t, x_v) SCRTBP(t, x_v, moon_p, C, main_p, sys_p), [0 TL1_ok/10], xx0L1_ok, options);

function U = U_ellips(x, y, z, sys_p)

a = sys_p(1);

b = sys_p(2);

c = sys_p(3);

G = sys_p(4);

rho = sys_p(5);

coeff1 = x^2 + y^2 + z^2;

coeff2 = x^2*(b^2+c^2) + y^2*(a^2+c^2) + z^2*(a^2+b^2);

coeff3 = (x^2*b^2*c^2 + y^2*a^2*c^2 + z^2*b^2*a^2 -1);

p = [coeff1 coeff2 coeff3];

small_k = max(roots(p));

w_small_k = asin(sqrt((a^2 - c^2) / (a^2 + small_k)));

k = (a^2 - b^2) / (a^2 - c^2);

F = ellipticF(w_small_k, k);

E = ellipticE(w_small_k, k);

fact1 = (2 * G * rho * pi * a * b * c) / sqrt(a^2 - c^2);

fact2 = 1 - x^2 / (a^2 - b^2) + y^2 / (a^2 - b^2);

fact3 = x^2 / (a^2 - b^2) - (a^2 - c^2) * y^2 / ((a^2 -b^2) * (b^2 - c^2)) + z^2 / (b^2 - c^2);

fact4 = (c^2 + small_k) * y^2 / (b^2 - c^2) - (b^2 + small_k * z^2 / (b^2 - c^2));

fact5 = sqrt(a^2 - c^2) / sqrt((a^2+ small_k) * (b^2 + small_k) * (c^2 + small_k));

U = fact1 * (fact2 * F + fact3 * E + fact4 * fact5);

end

function U = U_sph_harm(x, y, z, C, sys_p)

G = sys_p(1);

M = sys_p(2);

R = sys_p(3);

l = sys_p(4);

[lambda, phi, r] = cart2sph(x, y, z);

sum_U = 0;

C_count = 1;

for l_count = 1:l

for m_count = 0:l_count

l_count_st = num2str(2 * l_count);

m_count_st = num2str(2 * m_count);

lm_st = [l_count_st, m_count_st];

lm = str2double(lm_st);

Plm = ass_legendre(lm);

Plm = Plm(sin(phi));

sum_U = sum_U + (R/r)^(2*l) * C(C_count) * cos(2 * m_count * lambda) * Plm;

C_count = C_count + 1;

end

end

U = G * M / r * (1 + sum_U);

end

function U_scr3bp = gravity_pot_scr3bp(x, y, z, moon_p, C, main_p, sys_p)

OM = sys_p(1);

L = sys_p(2);

mu = sys_p(3);

U_scr3bp = (1/2) * (x^2 + y^2) + @(x, y, z) U_sph_harm((x + mu), y, z, C, main_p) + @(x, y, z) U_ellips((x + mu - 1), y, z, moon_p);

U_scr3bp = U_scr3bp / (OM^2 * L^2);

end

function dx_v = SCRTBP(~, x_v, moon_p, C, main_p, sys_p)

dx_v = zeros(6,1);

x = x_v(1); %rx

y = x_v(2); %ry

z = x_v(3); %rz

vx = x_v(4); %vx

vy = x_v(5); %vy

vz = x_v(6); %vz

syms x y z

U_scrbp_x = matlabFunction(diff(@(x, y, z) gravity_pot_scr3bp(x, y, z, moon_p, C, main_p, sys_p), x));

U_scrbp_y = matlabFunction(diff(@(x, y, z) gravity_pot_scr3bp(x, y, z, moon_p, C, main_p, sys_p), y));

U_scrbp_z = matlabFunction(diff(@(x, y, z) gravity_pot_scr3bp(x, y, z, moon_p, C, main_p, sys_p), z));

dx_v(1) = vx; %dx

dx_v(2) = vy; %dy

dx_v(3) = vz; %dz

dx_v(4) = U_scrbp_x + 2 * vy; %ddx

dx_v(5) = U_scrbp_y - 2 * vx; %ddy

dx_v(6) = U_scrbp_z; %ddz

end

function Plm = ass_legendre(lm)

switch lm

case 20

Plm = @(x) (1/2) * (3 * x^2 - 1);

case 22

Plm = @(x) 3 * (1 - x^2);

case 40

Plm = @(x) (1/8) * (35 * x^4 - 30 * x^2 + 3);

case 42

Plm = @(x) (15/2) * (7 * x^2 - 1) * (1 - x^2);

case 44

Plm = @(x) 105 * (1 - x^2);

end

[Dyuman Joshi]

Use end to close all the function definitions.

You have not defined the function/variable - ass_legendre. Thus we can not run your code to reproduce the error.

[Astrid Barzaghi]

My bad! edited to include ass_legendre

[Sam Chak]

Since the potential energy function U(x) is known, can you obtain the force at any position by taking the derivative of the potential using pen and paper? It feels solid if you can 'see' the equation rather than hiding it inside the matlabFunction().

[Debadipto]

Hi, what is the intended content of dx_v(4), dx_v(5) and dx_v(6)? A symbolic expression, function handle or a real value?