Cesario's answer is almost there, but I'd like to add a couple points of correction/improvement.

1) The problem formulation has constraints on the control (`0<=u<=1`) and positivity constraints on the states.  I would suggest rather than solving for root where the derivative is zero (equ3 above), use `Maximize` and add the constraints there:

```
parms = {s -> 10, m1 -> 0.02, m2 -> 0.5, m3 -> 4.4, r -> 0.03, 
   Tm -> 1500, k -> 0.000024, M -> 300, A -> 1, T0 -> 800, 
   Ti0 -> 0.04, V0 -> 1.5};
x1 = T[t];
x2 = Ti[t];
x3 = V[t];
X = {x1, x2, x3};
Lambda = {la1[t], la2[t], la3[t]};
U = u[t];
f1 = s/(1 + x3) - m1 x1 + r x1 (1 - (x1 + x2)/Tm) - U k x1 x3;
f2 = U k x3 x1 - m2 x2;
f3 = M m2 x2 - m3 x3;
F = {f1, f2, f3};
L = A x1 - (1 - U)^2;
H = L + Lambda . F;
equ1 = Thread[D[X, t] == Grad[H, Lambda]];
equ2 = Thread[D[Lambda, t] == -Grad[H, X]];
(*equ3=D[H,U]==0;
solU=Solve[equ3,U][[1]];*)
cons = {0 <= u[t] <= 1 && T[t] > 0 && Ti[t] > 0 && V[t] > 0, k > 0};
equ3 = Last@Maximize[{H, cons}, U];
solU = FullSimplify[equ3, cons]
```

This (nearly) yields the piecewise function reported in the original work.  Please note that the third ODE in the SE problem is missing a term found in the original problem, but the idea is the same.

[![ustar][1]][1]

2) The boundary conditions for lambda are given as lambda[0] in Cesario's answer, but these should apply to the end time tf, which is the gradient of the end cost (Mayer term) and here {0,0,0}.  Programmatically:

```
M=0; (*Mayer cost*)
tf=100;  (* from paper *)
cinitsLambda = Thread[(Lambda /. {t -> tf}) == Grad[M, X /. t -> tf]];
```
[![lambda terminal boundary condtions][2]][2]


Now the setup for NDSolve can be completed:
```
cinits = Join[cinitsX, cinitsLambda];
vars = Join[X, Lambda];
DE = Join[equs, cinits] /. solU /. parms;
```

3. And you might hope that the following NDSolve statement would work:

```
solDE = NDSolve[DE, vars, {t, 0, 0.01}];
Plot[Evaluate[vars /. solDE], {t, 0, 0.01}]
```

But alas it does not seem to converge in an overnight run, nor does the original equations from the paper. (vide supra)

I know that BVPs are tough problems to solve, so I tried an [alternative implementation][3] in my personal fork of PSOPT pseudospectral code, as well as the trapezoidal method in Mathematica.  These are *direct* solvers (vs. the Pontryagin indirect method) and also have problems with these equations. I have come to conclusion that there's something fishy or extreme stiffness.

I have tried this code on numerous multivariate optimal control problems and it seems to yield answers from Ross, Bryson and Ho, et al.


  [1]: https://i.sstatic.net/laKQr.png
  [2]: https://i.sstatic.net/Exn9y.png
  [3]: https://github.com/ecbrown/psopt/tree/examples/hiv