function optimal2(T1,T2)

% the problem is to find the final time for the system of the form:
% dX/dt = AX + Bu ; X(0) given
% such that J = 1/2X'(T)S(T)X(T) + 1/2int(rou + X'QX + ru^2, 0, T)
% be minimum
%
% this function plots the boundary condition difference for a range of
% final times between user defined values for T1 and T2. Then the answer
% can be found from the graph i.e. where the curve crosses the horizontal axis
%
% Ref.: Optimal Control by F. Leuis & V. Syrmous, 1995

% set the problem parameters

A=[0 1;0 0 ];
B=[0;1];
R=1;
Q=20*eye(2);
rou=4;
X0=[10;10];
k=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for T=T1:.01:T2
    
   k=k+1;
   [t,s]=ode45(@optimal2_ode,[T 0],[5 0 0 5]);

   s1=s(:,1);
   s2=s(:,2);
   s3=s(:,3);
   s4=s(:,4);

   ds1_t0=(s1(end-1)-s1(end))/(t(end-1)-t(end));
   ds2_t0=(s2(end-1)-s2(end))/(t(end-1)-t(end));
   ds3_t0=(s3(end-1)-s3(end))/(t(end-1)-t(end));
   ds4_t0=(s4(end-1)-s4(end))/(t(end-1)-t(end));

   ds_t0=[ds1_t0 ds2_t0; ds3_t0 ds4_t0];
   boundry_conditition_difference(k)=rou-X0'*ds_t0*X0;
   
end

plot([T1:.01:T2],boundry_conditition_difference)

function ds = optimal2_ode(t,s)

ds=zeros(4,1);

ds(1)=10-s(2)*s(3);
ds(2)=s(1)+10-s(2)*s(4);
ds(3)=s(1)+10-s(4)*s(3);
ds(4)=s(2)+s(3)+10-s(4)^2;