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];
    cinitsX = {x1 == T0, x2 == Ti0, x3 == V0} /. {t -> 0};
    cinitsLambda = Thread[Lambda == 0] /. {t -> 0};
    U = u[t];
    l1 = lambda1[t];
    l2 = lambda2[t];
    l3 = lambda3[t];
    f1 = s/(1 + x3) - m1 x1 + r x1 (1 - (x1 x2)/Tm);
    f2 = U k x3 x1 - m2 x2;
    f3 = M m2 x2 - m3 x3;
    X = {x1, x2, x3};
    F = {f1, f2, f3};
    Lambda = {l1, l2, l3};
    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]];
    equs = Join[equ1, equ2]
    cinits = Join[cinitsX, cinitsLambda]
    vars = Join[X, Lambda];
    DE = Join[equs, cinits] /. solU /. parms
    solDE = NDSolve[DE, vars, {t, 0, 1}]