how to figure out an Optimal control problem while occurring some system sensitivity?

I am solving a optimal control problem by using Gradient Decent Method in MAT LAB. The problem is while I put some small time it is giving me some output but while I increased the time the program runs but the Mat lab has encountered a crash report and it closed.

I want to know is the system is so sensitive or there are some bugs in my Code?

if true
% code 
function temon
clc
clear all
close all
alpha1=2;
beta1=-1.5;
gamma1=1;
alpha2=3.5;
beta2=-7;
gamma2=100;
eps = 1e-6;
step = 0.001;
time_factor=1;
t0 = 0; tf = 0.5;
t_segment = 20;

INPUT_VECTOR=[alpha1 beta1 gamma1 alpha2 beta2 gamma2 time_factor]; %EG2OC_Descent Example 2 of optimal control tutorial. % This example is from D.E.Kirk's Optimal control theory: an introduction % example 6.2-2 on page 338-339. % Steepest descent method is used to find the solution.

options = odeset('RelTol', 1e-4, 'AbsTol',[1e-4 1e-4]); % R = 0.1;

Tu = linspace(t0, tf, t_segment); % discretize time u = ones(1,t_segment); % guessed initial control u=1 initx = [0.5 6]; % initial values for states %initp = [0 0.001]; % initial values for costates max_iteration = 100; % Maximum number of iterations J_old=0 for i = 1:max_iteration % 1) start with assumed control u and move forward

u( u>10 )=10; u( u<0 )=0; [Tx,X] = ode45(@(t,x) stateEq(t,x,u,Tu,INPUT_VECTOR), [t0 tf], initx, options);

   % 2) Move backward to get the trajectory of costates
   x1 = X(:,1); x2 = X(:,2);

   p2T=2*gamma2*(X(end,2)-0);
   initp = [0 p2T];

   [Tp,P] = ode45(@(t,p) costateEq(t,p,u,Tu,x1,x2,Tx,INPUT_VECTOR), [tf t0], initp, options);
   p1 = P(:,1);
   p2 = P(:,2);

   % Important: costate is stored in reverse order. The dimension of
   % costates may also different from dimension of states
   % Use interploate to make sure x and p is aligned along the time axis
   p1 = interp1(Tp,p1,Tx);

   % Calculate gammaH with x1(t), x2(t), p1(t), p2(t)
   dH = pH(x1,p1,Tx,u,Tu,INPUT_VECTOR);
   H_Norm = dH'*dH;

   % Calculate the cost function

   J(i,1) = gamma2*(((X(end,2)^2)-0))

   %J(i,1) = tf*(((x1')*x1 + (x2')*x2)/length(Tx) + 0.1*(u*u')/length(Tu));

   % if dH/du < epslon, exit
   J_n = J_old - J(i,1)
   J_abs =  abs(J_n)
   % if dH/du < epslon, exit
   if H_Norm < eps
       % Display final cost
       J(i,1)
       break;
   elseif J_abs< eps
          J(i,1)
       break;
   else
       % adjust control for next iteration
       u_old = u;
       u = AdjControl(dH,Tx,u_old,Tu,step)
   end; 
   J_old=J(end,1);
end

% plot the state variables & cost for each iteration figure(1); grid on hold on %x1 = Tx(:,1); x2 = X(:,2); plot(Tx(:,1), x1 ,'g-'); plot(Tx(:,1), x2 ,'b-'); hold on; plot(Tu,u,'r'); legend('x1 = x1(t)','x2 = x2(t)','u*=u(t)'); xlabel('time','fontsize',16); ylabel('states','fontsize',16); set(gca,'FontSize',16); %ylim([-1 11]); hold off; print -djpeg90 -r300 eg2_descent.jpg

figure(2) grid on hold on plot(Tp(:,1), P(:,1) ,'r-'); plot(Tp(:,1), P(:,2) ,'b-'); legend('p1 = p1(t)','p2 = p2(t)'); xlabel('time','fontsize',16); ylabel('costates','fontsize',16); set(gca,'FontSize',18); hold off; print -djpeg90 -r300 eg2_descent.jpg

figure(3); grid on hold on plot(J,'r'); xlabel('Iteration number','fontsize',16); ylabel('J(u)= gamma2*(x2(T)-xd)^2','fontsize',16); legend('J(u)= gamma2*(x2(T)-xd)^2 '); f=J(1,1); st=(f/1.1);

s = strcat('Final cost is: J=',num2str(J(end,1)));

r= strcat('Final cost is: J=',num2str(J_abs)); text(1,f,s,'fontsize',16); text(st,f/1.1,r,'fontsize',16); set(gca,'FontSize',16); %xlim([-1 100000]); ylim([st f]); print -djpeg90 -r300 eg2_iteration.jpg

if i == max_iteration disp('Stopped before required residual is obtained.'); end

% State equations function dx = stateEq(t,x,u,Tu, INPUT_VECTOR)

alpha1 = INPUT_VECTOR(1,1); beta1 = INPUT_VECTOR(1,2); gamma1 = INPUT_VECTOR(1,3); alpha2 = INPUT_VECTOR(1,4); beta2 = INPUT_VECTOR(1,5); gamma2 = INPUT_VECTOR(1,6); time_factor = INPUT_VECTOR(1,7);

dx = zeros(2,1); u = interp1(Tu,u,t); % Interploate the control at time t %dx(1) = -2*(x(1) + 0.25) + (x(2) + 0.5)*exp(25*x(1)/(x(1)+2)) - (x(1) + 0.25).*u; %dx(2) = 0.5 - x(2) -(x(2) + 0.5)*exp(25*x(1)/(x(1)+2)); dx(1) = time_factor*(alpha1*x(1) + beta1*x(1)*x(2) + (gamma1*x(1)).*u); dx(2) = time_factor*(alpha2*x(2) + beta2*x(1)*x(2));

% Costate equations function dp = costateEq(t,p,u,Tu,x1,x2,xt,INPUT_VECTOR)

alpha1 = INPUT_VECTOR(1,1); beta1 = INPUT_VECTOR(1,2); gamma1 = INPUT_VECTOR(1,3); alpha2 = INPUT_VECTOR(1,4); beta2 = INPUT_VECTOR(1,5); gamma2 = INPUT_VECTOR(1,6); time_factor = INPUT_VECTOR(1,7);

dp = zeros(2,1); x1 = interp1(xt,x1,t); % Interploate the state varialbes x2 = interp1(xt,x2,t); u = interp1(Tu,u,t); % Interploate the control

dp(1)= time_factor*(-(p(1).*(alpha1 + beta1*x2 + gamma1.*u) + p(2).*(beta2*x2))); dp(2)= time_factor*(-(p(1).*(beta1*x1) + p(2).*(alpha2*+ beta2*x1)));

% Partial derivative of H with respect to u function dH = pH(x1,p1,tx,u,Tu, INPUT_VECTOR) alpha1 = INPUT_VECTOR(1,1); beta1 = INPUT_VECTOR(1,2); gamma1 = INPUT_VECTOR(1,3); alpha2 = INPUT_VECTOR(1,4); beta2 = INPUT_VECTOR(1,5); gamma2 = INPUT_VECTOR(1,6); time_factor = INPUT_VECTOR(1,7);

% interploate the control u = interp1(Tu,u,tx); dH = p1.*(x1 * gamma1*time_factor);

% Adjust the control function u_new = AdjControl(pH,tx,u,tu,step)

% interploate dH/du pH = interp1(tx,pH,tu); u_new = u - step*pH; u_new (u_new >10 )=10; u_new (u_new <0 )=0;

end