Computing at each time step a component of a predefined vector in a function called by an ode45

Hi everyone!

I wonder if anyone knows how to solve this issue:

I want to compute in an embedded ode function, components from vectors which I get from .txt files. I have tried it using a loop within the function, from 1 to the length of the vectors I want to include (all vectors have the same length). As the ode45 solver, or ode113 use an adaptive time step, I don't know if this method is correct.

I attach the code, in order to help you to understand the problem.

Thank you in advance

%% Introduccion de datos necesarios para la perturbacion

% DATOS PERTURBACION J2
J2 = 2e-7;
Rs = 696000;
mu1 = 132.712e9; % parametro gravitacional Sol (km^3/s^2)

% DATOS PERTURBACION PRESION DE RADIACION SOLAR
c = 2.99e8; % velocidad de la luz (m/s)
% S = 1367;   % constante solar (w/m^2)
SO = 63.15e6; % Intensidad de potencia de radiacion (W/m^2)
As = 4.474;  % Area frontal del escudo termico (m^2)
mpsp = 605; % Masa en seco + payload (kg)
CR = 2;     % Coeficiente de presion de radiacion (1 < CR < 2)
nu = 1;


% DATOS PERTURBACION SOL 

% Para la prueba con coordenadas baricentricas

%rSolx = CoordSol(:,4);
%rSoly = CoordSol(:,5);
%rSolz = CoordSol(:,6);

This are the vectors I get from .txt files.

% DATOS PERTURBACION MERCURIO
mu_m = 22030;
rMx = CoordMerc(:,4);
rMy = CoordMerc(:,5);
rMz = CoordMerc(:,6);

% DATOS PERTURBACION VENUS
mu_v = 324900; % parametro gravitacional Venus (km^3/s^2)
rVx = CoordVenus(:,4);
rVy = CoordVenus(:,5);
rVz = CoordVenus(:,6);

% DATOS PERTURBACION TIERRA
mu_t = 398600;
rTx = CoordTierra(:,4);
rTy = CoordTierra(:,5);
rTz = CoordTierra(:,6);

% DATOS PERTURBACION MARTE
mu_mr = 42828;
rMrx = CoordMarte(:,4);
rMry = CoordMarte(:,5);
rMrz = CoordMarte(:,6);

% DATOS PERTURBACION JUPITER
mu_j = 126686000;
rJx = CoordJupiter(:,4);
rJy = CoordJupiter(:,5);
rJz = CoordJupiter(:,6);

% DATOS PERTURBACION SATURNO
mu_s = 37931000;
rSx = CoordSaturno(:,4);
rSy = CoordSaturno(:,5);
rSz = CoordSaturno(:,6);

% DATOS PERTURBACION URANO
mu_u = 5794000;
rUx = CoordUrano(:,4);
rUy = CoordUrano(:,5);
rUz = CoordUrano(:,6);

% DATOS PERTURBACION NEPTUNO
mu_n = 6835100;
rNx = CoordNeptuno(:,4);
rNy = CoordNeptuno(:,5);
rNz = CoordNeptuno(:,6);


%% METODOS DE INTEGRACION

% Vector estado
y0 = [r0 v0]';

% =========================================
% =========================================

% Paso temporal. Debe ser en segundos, en funcion de la precision que se
% desee utilizar, el paso temporal sera:

% h = 1; Precision de 1 segundo
% h = 60; Precision de 1 minuto
% h = 3600; Precision de 1 hora
% h = 86400; Precision de 1 dia


h = 86400;
tspan = t0:h:tf;

% Metodo de paso adaptativo RK45

% [t,y1] = rkf45(@rates,tspan,y0,1e-10);
 
% Metodo ode45
 [t,y1] = ode45(@rates, tspan, y0,1e-10);

% Metodo ode113
% [t,y1] = ode113(@rates, tspan, y0,1e-10);


% =========================================
% =========================================

% Metodo de paso constante
% rk = 4;
% [t, y1] = rk1_4(@rates, [t0 tf], y0, h, rk);

% =========================================
% =========================================
% Predictor - corrector
%[t,y1] = heun(@rates, [t0 tf], y0, h);

% =========================================
% =========================================


%% CALCULO DE LA ACELERACION


function [ dydt ] = rates(t,f)
%{
Esta función se encarga de calcular la aceleración del vector 
t - tiempo
f - vector columna que contiene el vector posición y el vector
    velocidad para cada tiempo t
x, y, z - componentes del vector posición r
r - módulo del vector posición r
vx, vy, vz - componentes del vector velocidad v
ax, ay, az - componentes del vector aceleración a
dydt - vector columna que contiene las componentes de los vectores de
velocidad y aceleración
%}
for i=1:length(rVx) 
    
x = f(1);
y = f(2);
z = f(3);
vx = f(4);
vy = f(5);
vz = f(6);
r = norm([x y z]);

% Calculo vector posicion relativa de Mercurio al satelite
rxm_s = rMx(i) - x;
rym_s = rMy(i) - y;
rzm_s = rMz(i) - z;
rM_s = norm([rxm_s rym_s rzm_s ]);
rM = norm([rMx(i)  rMy(i) rMz(i) ]);


% Calculo vector posicion relativa de Venus al satelite
rxv_s = rVx(i) - x;
ryv_s = rVy(i) - y;
rzv_s = rVz(i) - z;
rV_s = norm([rxv_s ryv_s rzv_s ]);
rV = norm([rVx(i)  rVy(i) rVz(i) ]);

% Calculo vector posicion relativa de La Tierra al satelite
rxt_s = rTx(i) - x;
ryt_s = rTy(i) - y;
rzt_s = rTz(i) - z;
rT_s = norm([rxt_s ryt_s rzt_s ]);
rT = norm([rTx(i)  rTy(i) rTz(i) ]);

% Calculo vector posicion relativa de Marte al satelite
rxmr_s = rMrx(i) - x;
rymr_s = rMry(i) - y;
rzmr_s = rMrz(i) - z;
rMr_s = norm([rxmr_s rymr_s rzmr_s ]);
rMr = norm([rMrx(i)  rMry(i) rMrz(i) ]);

% Calculo vector posicion relativa de Saturno al satelite
rxs_s = rSx(i) - x;
rys_s = rSy(i) - y;
rzs_s = rSz(i) - z;
rS_s = norm([rxs_s rys_s rzs_s ]);
rS = norm([rSx(i)  rSy(i) rSz(i) ]);

% Calculo vector posicion relativa de Jupiter al satelite
rxj_s = rJx(i) - x;
ryj_s = rJy(i) - y;
rzj_s = rJz(i) - z;
rJ_s = norm([rxj_s ryj_s rzj_s ]);
rJ = norm([rJx(i)  rJy(i) rJz(i) ]);

% Calculo vector posicion relativa de Urano al satelite
rxu_s = rUx(i) - x;
ryu_s = rUy(i) - y;
rzu_s = rUz(i) - z;
rU_s = norm([rxu_s ryu_s rzu_s ]);
rU = norm([rUx(i)  rUy(i) rUz(i) ]);

% Calculo vector posicion relativa de Neptuno al satelite
rxn_s = rNx(i) - x;
ryn_s = rNy(i) - y;
rzn_s = rNz(i) - z;
rN_s = norm([rxn_s ryn_s rzn_s ]);
rN = norm([rNx(i)  rNy(i) rNz(i) ]);

% norma

% rSol = norm([rSolx(i) rSoly(i) rSolz(i)]); % Caso de baricentricas
 
% Calculo aceleraciones con perturbaciones
ax = -mu*(x)/(r)^3 + ((3/2)*(J2*mu1*Rs^2/r^4)*((x/r)*(5*(z^2/r^2)-1))) ...
     +(-nu*((SO*(Rs/(r))^2)/c)*(CR*As/mpsp))*(x)/(r)^3 + ...
     mu_m*(rxm_s/rM_s^3-rMx(i)/rM^3) + ...
     mu_v*(rxv_s/rV_s^3-rVx(i)/rV^3) +...
     mu_t*(rxt_s/rT_s^3-rTx(i)/rT^3)+...
     mu_mr*(rxmr_s/rMr_s^3-rMrx(i)/rMr^3) + ...
     mu_j*(rxj_s/rJ_s^3-rJx(i)/rJ^3) + ...
     mu_s*(rxs_s/rS_s^3-rSx(i)/rS^3) + ...
     mu_u*(rxu_s/rU_s^3-rUx(i)/rU^3) + ...
     mu_n*(rxn_s/rN_s^3-rNx(i)/rN^3);
       
ay = -mu*(y)/(r)^3 + ((3/2)*(J2*mu1*Rs^2/r^4)*((y/r)*(5*(z^2/r^2)-1))) ...
    + (-nu*((SO*(Rs/(r))^2)/c)*(CR*As/mpsp))*(y)/(r)^3 + ... 
     mu_m*(rym_s/rM_s^3-rMy(i)/rM^3) + ...
     mu_v*(ryv_s/rV_s^3-rVy(i)/rV^3) +...
     mu_t*(ryt_s/rT_s^3-rTy(i)/rT^3)+...
     mu_mr*(rymr_s/rMr_s^3-rMry(i)/rMr^3) + ...
     mu_j*(ryj_s/rJ_s^3-rJy(i)/rJ^3) + ...
     mu_s*(rys_s/rS_s^3-rSy(i)/rS^3) + ...
     mu_u*(ryu_s/rU_s^3-rUy(i)/rU^3) + ...
     mu_n*(ryn_s/rN_s^3-rNy(i)/rN^3);
        
az = -mu*(z)/(r)^3 + ((3/2)*(J2*mu1*Rs^2/r^4)*((z/r)*(5*(z^2/r^2)-3))) ...
    + (-nu*((SO*(Rs/(r))^2)/c)*(CR*As/mpsp))*(z)/(r)^3 + ...
     mu_m*(rzm_s/rM_s^3-rMz(i)/rM^3) + ...
     mu_v*(rzv_s/rV_s^3-rVz(i)/rV^3) +...
     mu_t*(rzt_s/rT_s^3-rTz(i)/rT^3)+...
     mu_mr*(rzmr_s/rMr_s^3-rMrz(i)/rMr^3) + ...
     mu_j*(rzj_s/rJ_s^3-rJz(i)/rJ^3) + ...
     mu_s*(rzs_s/rS_s^3-rSz(i)/rS^3) + ...
     mu_u*(rzu_s/rU_s^3-rUz(i)/rU^3) + ...
     mu_n*(rzn_s/rN_s^3-rNz(i)/rN^3);

dydt = [vx vy vz ax ay az]';

end

end
end