Solving Differential Equation System for HIV Treatment Model

I was working on a project about optimal strategies for HIV treatment, models used from [Butler, Kirschner, and Lenhart] 1997. This model explains the spread of HIV viruses in the human body, where there is one control function u(t). My work is following pontryagin maximum principle. But i have a problem solving the differential equation system, where there are 6 differential equations with 6 initial conditions. Here i use software Mathematica 11 when I execute the code there is no error that appears and he just returning DSolve code that i write.

ClearAll["Global`*"];
dT = s/(1 + V[t]) - m1 T[t] + r T[t] (1 - (T[t] + Ti[t])/Tm) -
u[t] k V[t] T[t];
dTi = u[t] k V[t] T[t] - m2 Ti[t];
dV = M m2 Ti[t] - m3 V[t];
H = A T[t] - (1 - u[t])^2 + l1[t] dT + l2[t] dTi + l3[t] dV ;
u[t] = 1/2 (2 - k l1[t] T[t] V[t] + k l2[t] T[t] V[t])
eq1 = D[l1[t], t] == D[H, l1[t]];
eq2 = D[l2[t], t] == D[H, l2[t]];
eq3 = D[l3[t], t] == D[H, l3[t]];
eq4 = D[T[t], t] == -D[H, T[t]];
eq5 = D[Ti[t], t] == -D[H, Ti[t]];
eq6 = D[V[t], t] == -D[H, V[t]];
s = 10; m1 = 0.02; m2 = 0.5; m3 = 4.4; r = 0.03; Tm = 1500; k = \
0.000024; M = 300; A = 1;

DSolve[{eq1, eq2, eq3, eq4, eq5, eq6, T == 800, Ti == 0.04,
V == 1.5, l1 == 0, l2 == 0, l3 == 0}, {l1[t], l2[t],
l3[t], T[t], Ti[t], V[t]}, t]

I cannot guarantee that the analytical solution exists, but is there something wrong with code I wrote? Or is there an alternative solution to complete the differential equation system?

I hope you are pleased to check the file that I attached. Control-Optimal.nb

Many thanks • Thanks for your comment, it give this error: NDSolve::deqn: Equation or list of equations expected instead of eq1 in the first argument {eq1,eq2,eq3,eq4,eq5,eq6,T==800,Ti==0.04,V==1.5,l1==0,l2==0,l3==0}. – Wahyu Surya Ningrat Dec 24 '18 at 13:25
• I really appreciate it, thank you – Wahyu Surya Ningrat Dec 24 '18 at 13:32
• I'm sure there is no analytical solution, so you'll need to use NDSolve. I noticed that the equations for T'[t], Ti'[t] and V'[t] in your code do not match the system you want to solve. In any case, you might want to include all your code in the question in case people do not want to download a notebook. – Chris K Dec 24 '18 at 14:42
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Dec 24 '18 at 14:46
• Oh i see, i will focus on NDSolve, thanks for your advice @ChrisK – Wahyu Surya Ningrat Dec 24 '18 at 23:27

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][];
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, 1}];
Plot[Evaluate[X /. solDE], {t, 0, 1}]
• thankyou very much, its help me to learn this algorithm a lot – Wahyu Surya Ningrat Jan 3 at 17:52