ODE45 with fmincon optimization error

options = optimoptions('fmincon','Display','iter','Algorithm','sqp');

x0=[400 31.095 26.718 19.356 13.096 6.49 0.758 -10.614 -27.413 -27.413];

A = [];
b = [];
Aeq = [];
beq = []; 
lb = zeros(1,10);
ub = zeros(1,10);
lb(1,1)=0;
ub(1,1)=10000;
for i=2:10
    lb(1,i)=-90;
    ub(1,i)=90;
end
lb
ub
global m_pay
global mis
mis=400*1000;
global phase_data
phase_data=[0 10 20 30 50 75 100 140 175 191.34 193 195 220 260 310 350 390 420 490 560 605;
    90 90 89.973 89.898 84.76 71.349 56.021 38.934 30.622 37.74 37.74 37.74 31.095 26.718 19.356 13.096 6.49 0.758 -10.614 -27.413 -27.413];

global m_s
m_s=[2717 544];
global m_p
m_p=[39463 4037];
global I_sp
I_sp=[304 325];
global m_dot
global T
global A_e
global r_e
global mu
global S_ref
global tout yout y0
%m_pay=400; %mass of payload
T=1000*[62.5,3.2]; %thrust of rocket
m_dot=T./I_sp; %thrust differential
T=9.806*T;
A_e=[0.9 0.2]; %area of rocket
r_e=6378000; %Earth radius
S_ref=[3.14 1]; %reference area
mu=398600*10^9; %mu
global pd
pd=[0:2000:86000;
    1.013e+5 7.950e+4 6.166e+4 4.722e+4 3.565e+4 2.650e+4 1.940e+4 1.417e+4 1.035e+4 7.565e+3 5.529e+3 4.047e+3 2.972e+3 2.188e+3 1.616e+3 1.197e+3 8.890e+2 6.634e+2 4.985e+2 3.771e+2 2.871e+2 2.200e+2 1.695e+2 1.313e+2 1.023e+2 7.977e+1 6.221e+1 4.833e+1 3.736e+1 2.872e+1 2.196e+1 1.669e+1 1.260e+1 9.459e+0 7.051e+0 5.220e+0 3.835e+0 2.800e+0 2.033e+0 1.467e+0 1.052e+0 7.498e-1 5.308e-1 3.732e-1;
    1.225e+0 1.007e+0 8.193e-1 6.601e-1 5.258e-1 4.135e-1 3.119e-1 2.279e-1 1.665e-1 1.216e-1 8.891e-2 6.451e-2 4.694e-2 3.426e-2 2.508e-2 1.841e-2 1.355e-2 9.887e-3 7.257e-3 5.366e-3 3.995e-3 2.995e-3 2.259e-3 1.714e-3 1.317e-3 1.027e-3 8.055e-4 6.389e-4 5.044e-4 3.962e-4 3.096e-4 2.407e-4 1.860e-4 1.429e-4 1.091e-4 8.281e-5 6.236e-5 4.637e-5 3.430e-5 2.523e-5 1.845e-5 1.341e-5 9.690e-6 6.955e-6
    340.3 332.5 324.6 316.5 308.1 299.5 295.1 295.1 295.1 295.1 295.1 296.4 297.7 299.1 300.4 301.7 303.0 306.5 310.1 313.7 317.2 320.7 324.1 327.5 329.8 329.8 328.8 325.4 322.0 318.6 315.1 311.5 308.0 304.4 300.7 297.1 293.4 290.7 288.0 285.3 282.5 279.7 276.9 274.1];
%first line altitude, second line pressure, third line density
global mach
mach=[0,0.7,0.95,1.1,1.2,1.4,1.7,2.5,4,50;
    0.35,0.45,0.75,1.1,1.15,1.05,0.9,0.65,0.55,0.55];
%first line mach number, second line drag coefficient
tout=[];
yout=[];
y0=[0 r_e 0 0 m_pay+sum(m_p)+sum(m_s)];
global iPhase;

for iPhase=1:12
    t0 = phase_data(1,iPhase);
    tf = phase_data(1,iPhase+1);
    [t, y] = ode45(@dynamics, [t0 tf], y0);
    tout = [tout; t];
    yout = [yout; y];
    plot (tout,yout)
    y0 = y(end,:);
    if iPhase == 10
        y0(5)=m_s(1,2)+m_p(1,2)+m_pay;
    end
    
    ss=yout(:,1:2);
    als=ss.*ss;
    al=als(:,1)+als(:,2);
    al=sqrt(al);
    plot(tout,(al-r_e)/1000)
end


[x2, fval2] = fmincon(@part2_1, x0, A, b, Aeq, beq, lb, ub, @part2_2, options);
%fun=@(x) -x(1);
%[x2, fval2] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, @part2_2, options);
x2
function [c,ceq]=part2_2(xx)
    global phase_data
    global m_pay m_s m_p
    global iPhase
    global tout yout
    global y0
    global r_e
    global mis mu
    phase_data(2,13:21)=xx(1,2:10);
    m_pay=xx(1,1);
    for iPhase=13:20
        t0 = phase_data(1,iPhase);
        tf = phase_data(1,iPhase+1);
        [t, y] = ode45(@dynamics, [t0 tf], y0);
        tout = [tout; t];
        yout = [yout; y];
        y0 = y(end,:);
        if iPhase == 10
            y0(5)=m_s(1,2)+m_p(1,2)+m_pay;
        end
    end
    r=sqrt(sum(y0(1:2).*y0(1:2))); %|r|
    v=sqrt(sum(y0(3:4).*y0(3:4))); %|v|
    ceq(1)=(r_e+mis-r);
    ceq(2)=(v-sqrt(mu/(r_e+mis)));
    ceq(3)=(sum(y0(1:2).*y0(3:4)));
    c=[];
end
%{
ss=yout(:,1:2);
als=ss.*ss;
al=als(:,1)+als(:,2);
al=sqrt(al);
plot(tout,(al-r_e)/1000)
%}


function out1 = part2_1(xx)

    
    x1=xx(1);
    out1=-x1;

end



function y = dynamics(t,x)
    global iPhase
    global m_pay
    global m_s
    global m_p
    global m_dot
    global phase_data
    global r_e
    global mach
    global S_ref
    global pd
    global T
    global A_e
    global mu
    r=sqrt(sum(x(1:2).*x(1:2))); %|r|
    v=sqrt(sum(x(3:4).*x(3:4))); %|v|
    alt=r-r_e; %altitude
    p=interp1(pd(1,:),pd(2,:),alt,'pchip');
    d=interp1(pd(1,:),pd(3,:),alt,'pchip');
    a=interp1(pd(1,:),pd(4,:),alt,'pchip');
    beta=looking(phase_data(1,:),phase_data(2,:),t);
    beta=beta*pi/180;
    y=[0 0 0 0 0];
    if iPhase<=9
        m=m_pay+sum(m_s)+sum(m_p)-m_dot(1,1)*t;
        y(5)=-m_dot(1,1);
        tr0=T(1,1)-A_e(1,1)*p;
        tr=tr0*[cos(beta) sin(beta)];
        C_d=interp1(mach(1,:),mach(2,:),v/a,'pchip');
        s=S_ref(1);
        
    end
    if iPhase==10
        m=m_pay+sum(m_s)+sum(m_p)-m_dot(1,1)*phase_data(1,10);
        y(5)=0;
        tr0=0;
        tr=[0 0];
        C_d=interp1(mach(1,:),mach(2,:),v/a,'pchip');
        s=S_ref(1);
        
    end
    if iPhase==11
        m=m_s(1,2)+m_p(1,2)+m_pay;
        y(5)=0;
        s=S_ref(2);
        tr0=0;
        tr=[0 0];
        C_d=0;
        
    end
    if iPhase>=12
        m=m_pay+m_s(1,2)+m_p(1,2)-m_dot(1,2)*(t-phase_data(1,12));
        y(5)=-m_dot(1,2);
        s=S_ref(2);
        tr0=T(1,2)-A_e(1,2)*p;
        tr=tr0*[cos(beta) sin(beta)];
        C_d=0;
        
    end
    g=-mu/(r.^3)*x(1:2);
    d=-0.5*d*v*s*C_d*x(3:4);
    f=g+(d+tr')/m;
    y(1:2)=x(3:4);
    y(3:4)=f;
    y=y';    
end

function res=looking(x,y,k)
    
    num = size(x);
    num=num(2)-1;
    jj=1;
    while x(1,jj)<=k
        jj=jj+1;
        if jj>num
            break
        end
    end
    jj=jj-1;
    res=y(jj)+(k-x(jj))/(x(jj+1)-x(jj))*(y(jj+1)-y(jj));
end

I am trying to optimize function with fmincon, but It keep shows warning and return initial value, not optimized value.

Error message also says 'Unable to meet integration tolerances without reducing the step size below the smallest value allowed'. What should I do to solve this error?