Error in ODE45 using fmincon, Arrays incompatible sizes

I am trying to integrate equations using ode45 while at the same time optimizing parameters within these. Input to fmincon are some parameters for functions within the equations to be integrated. In the objective_fun_circ function, the objective function is calculated based on the results of integration with ode45, which is done by launching a second function called trajectory_optimization.

% Initialization control variables
eta_min     = 0.01;                                       
eta_list0   = [0.99 0.99]      
function_list0 = [theta_01 theta_f1 theta_f2 a1 a2 csi1 csi2]; %values are defined previously to this vector I have not reported the
%entire code (same for data)
% Initial guess in vector form
x0 = [eta_list0(:);function_list0(:)]
% Lower and upper boundaries of control variables
function_lb = [-pi/2 -pi/2 -pi/2 60 60 -1 -1]';
function_ub = [pi/2 pi/2 pi/2 140 140 1 1]';
lb = [eta_min*ones(data.n_stage,1);function_lb];
ub = [ones(data.n_stage,1);function_ub];
%Optimization
options = optimoptions('fmincon','Display','iter-detailed', 'MaxFunctionEvaluations', 2000);
[x,fval,exitflag,output,lambda,grad,hessian] = fmincon(@(x) objective_fun_circ(x,data, Target),x0,[],[],[],[],lb,ub,@(x) mycon(x,data, Target),options)

In the objective_fun_circ function, the objective function is calculated based on the results of integration with ode45, which is done by launching a second function called trajectory_optimization.

%% Reshape the control variable
eta_list = reshape(x(1:data.n_stage),[data.n_stage, 1]);
function_list  = reshape(x(data.n_stage+1:end),[7,1]);
    
p0 = trajectory_optimized(eta_list,function_list,data);

ODE45 is finally used in this function.

function [p0,p_total,t_total] = trajectory_optimized(eta_list,function_list,data)
global g0_Earth g0_Mars mu_mars R_mars
p0 = [r_0; th_0; v_0; gamma_0; m_0; d_loss; g_loss; delta_v]; 
%values are defined previously to this vector I have not reported the
%entire code (same for "data")
t0=0;
t_total = [];       
p_total = [];      

for stage = 1:data.N_stage 
    p0(5) = data.m0(stage);
    stage_time = data.isp*g0_Earth*data.fuel/data.thrust/sum(eta_list(stage));
    [t,p] = ode45(@(t,p) EOM_optimized(t,p,eta_list(stage),function_list,stage,data),[0.01,eta_list(stage)*stage_time],p0,options);
        
       
        p_total = [p_total;p];          
        t_total = [t_total;t0+t];      
        
        t0 = t_total(end);
        p0 = p_total(end,:)';
    
end
end

In EOM_optimized equations to be integrated are written

function dp = EOM_optimized(t,p,eta,function_list,stage,data)
global g0_Earth g0_Mars mu_mars R_mars
thrust = data.thrust;
isp = data.isp;
Cd = data.cd;
Area = data.Area;
drag=data.drag
r = p(1);
theta = p(2);
v = p(3);
gamma = p(4);
m = p(5);
d_loss = p(6);
g_loss = p(7);
delta_v = p(8);

...
    
dr = v*sin(gamma);
dtheta = v*cos(gamma)/r;
if (r-R_mars) <= 250
dv = (eta*thrust/m)-0.5*rho/m*Area*Cd*v^2-(mu_mars/r^2)*sin(gamma);
else
    if stage==1
    dv = (eta*thrust/m)*cos(data.alpha(t,function_list)-gamma)-0.5*rho/m*Area*Cd*v^2-(mu_mars/r^2)*sin(gamma);
    else
    dv = (eta*thrust/m)*cos(data.alpha2(t,function_list)-gamma)-0.5*rho/m*Area*Cd*v^2-(mu_mars/r^2)*sin(gamma);
    end
end

if (r-R_mars) <= 250
dgamma=0;
else
    if stage==1
    dgamma=(eta*thrust/m)*sin(data.alpha(t,function_list)-gamma)-(1/v)*(mu_mars/r^2)*cos(gamma)+(v*cos(gamma))/r;
    else
    dgamma=(eta*thrust/m)*sin(data.alpha2((t),function_list)-gamma)-(1/v)*(mu_mars/r^2)*cos(gamma)+(v*cos(gamma))/r;
    end
end

dm = -eta*thrust/(g0_Earth*isp);
d_delta_v = isp*g0_Earth*log(m/(m+dm));
d_drag_loss = drag/m;
d_gravity_loss = (g0_Mars/(r/R_mars)^2)*sin(gamma);

dp = [dr;dtheta; dv;dgamma;dm;d_drag_loss; d_gravity_loss; d_delta_v];
end

An error is reported in ode45 regarding incompatible sizes for arrays in the calculation of “ynew”, the problem is that by running the code without fmincon or allowing only a few iterations (2-3) the error does not occur, so I cannot understand the problem. Also by changing some input variables in present in “data” or by changing the initial integration time of ode45 by very little (e.g. from 0 to 0.01) the error changes (e.g. the operation that fails involves y4 or y5 instead of ynew) or disappears, I however need to leave the current values in “data”.

[Torsten]

Make it easy for us and include executable code that reproduces the error.