% script_BrysonProblem_AnalyticSolution
% thanks to Martin Neuenhofen for providing script
t0      = 0;
tE      = 1.5;
t1      = 0.5/tanh(1.5);

% t2 is unique solution of 
% (0.5+t1-0.5*t1^2)*exp(t1-t2)==te-t2+0.5*(tE-t2)^2
% in interval [t0,tE].
t2      = 1.05462;

figure;
hold on;
box on;
axis manual;
axis([t0,tE,-1.5,1.5]);
grid on;

% interval [t0,t1]
u       = @(t)  1+0*t;
x2      = @(t)  -t;
x1      = @(t)  0.5+t-0.5*t.^2;

t       = linspace(t0,t1,1000);
plot(t,x1(t),'blue','linewidth',1);
plot(t,x2(t),'red','linewidth',1);
plot(t,u(t) ,'black','linewidth',1);

% interval [t1,t2]
x1_0    = x1(t1);
x2_0    = x2(t1);
x1      = @(t)  x1_0 ./ exp(t-t1);
x2      = @(t)  x1_0 * sinh(t-t1) + x2_0 * exp(t-t1);
u       = @(t)  -x1(t)-x2(t);

t       = linspace(t1,t2,1000);
plot(t,x1(t),'blue','linewidth',1);
plot(t,x2(t),'red','linewidth',1);
plot(t,u(t) ,'black','linewidth',1);

% interval [t2,tE]
u       = @(t)  -1+0*t;
x2      = @(t)  -1.5+t;
x1      = @(t)  (tE-t)+0.5*(tE-t).^2;

t       = linspace(t2,tE,1000);
plot(t,x1(t),'blue','linewidth',1);
plot(t,x2(t),'red','linewidth',1);
plot(t,u(t) ,'black','linewidth',1);

legend({'$x_1(t)$','$x_2(t)$','$u(t)$'},...
    'Interpreter','LateX','fontsize',12);
title('Bryson Singular Arc Optimal Control Solution',...
    'Interpreter','LateX','fontsize',12);
xlabel('time $t$','Interpreter','LateX','fontsize',12);
ylabel('states $x_1(t),x_2(t)$\,, control $u(t)$',...
    'Interpreter','LateX','fontsize',12);