Structuring the script. 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]]; equs = Join[equ1, equ2]; cinitsX = Thread[X == {T0, Ti0, V0}] /. {t -> 0}; cinitsLambda = Thread[Lambda == 0] /. {t -> 0}; cinits = Join[cinitsX, cinitsLambda]; vars = Join[X, Lambda]; DE = Join[equs, cinits] /. solU /. parms; solDE = NDSolve[DE, vars, {t, 0, 0.01}]; Plot[Evaluate[vars /. solDE], {t, 0, 0.01}]