with(VectorCalculus);
with(linalg);


pi := 4; eta := 6; mh := .1; mv := .3; bh := .15; Nh := 400; Nv := 200; bv := .25; b := .8; d := .1; p := .5; B := 10; A := 1; alpha := .25; beta := .7; K := 400; r := .5; c := .9; q := .8; Sh0 := 225; Sv0 := 100; Ih0 := 175; Iv0 := 600; P0 := 50; T := 35; B := 10;
                  


eq1 := diff(L1(t), t) = -L1(t)*(-bh*b*Iv(t)/Nh-mh)-L2(t)*bh*b*Iv(t)/Nh; eq2 := diff(L2(t), t) = -L2(t)*(-mh-d)+L3(t)*bv*b*Sv(t)/Nv-L4(t)*bv*b*Sv(t)/Nv; eq3 := diff(L3(t), t) = -(L4(t)*Iv(t)+L3(t)*Sv(t))*L3(t)/(2*B^2)-L3(t)*(-bv*b*Ih(t)/Nv-mv-L3(t)*Sv(t)/(2*B^2)-(L4(t)*Iv(t)+L3(t)*Sv(t))/(2*B^2)-P(t)*alpha)-L4(t)*(bv*b*Ih(t)/Nv-L3(t)*Iv(t)/(2*B^2))-L5(t)*P(t)*alpha; eq4 := diff(L4(t), t) = -A-(L4(t)*Iv(t)+L3(t)*Sv(t))*L4(t)/(2*B^2)+L1(t)*bh*b*Sh(t)/Nh-L2(t)*bh*b*Sh(t)/Nh+L3(t)*L4(t)*Sv(t)/(2*B^2)-L4(t)*(-mv-L4(t)*Iv(t)/(2*B^2)-(L4(t)*Iv(t)+L3(t)*Sv(t))/(2*B^2)-P(t)*alpha)-L5(t)*P(t)*alpha; eq5 := diff(L5(t), t) = L3(t)*Sv(t)*alpha+L4(t)*Iv(t)*alpha-L5(t)*(r*(1-P(t)/Nv)-P(t)*r/Nv+alpha*(Sv(t)+Iv(t))-q*c); eq6 := diff(Sh(t), t) = pi-bh*b*Sh(t)*Iv(t)/Nh-mh*Sh(t); eq7 := diff(Ih(t), t) = bh*b*Sh(t)*Iv(t)/Nh-mh*Ih(t)-d*Ih(t); eq8 := diff(Sv(t), t) = eta-bv*b*Sv(t)*Ih(t)/Nv-mv*Sv(t)-P(t)*Sv(t)*alpha-L3(t)*Sv(t)^2/(2*B^2)-L4(t)*Sv(t)*Iv(t)/(2*B^2); eq9 := diff(Iv(t), t) = -L3(t)*Sv(t)*Iv(t)/(2*B^2)+bv*b*Sv(t)*Ih(t)/Nv-mv*Iv(t)-P(t)*Iv(t)*alpha-L4(t)*Iv(t)^2/(2*B^2); eq10 := diff(P(t), t) = P(t)*r*(1-P(t)/Nh)+P(t)*alpha*(Sv(t)+Iv(t))-q*P(t)*c

fcns := {Ih(t), Iv(t), L1(t), L2(t), L3(t), L4(t), L5(t), P(t), Sh(t), Sv(t)}; a := dsolve({eq1, eq10, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, Ih(0) = Ih0, Iv(0) = Iv0, L1(T) = 0, L2(T) = 0, L3(T) = 0, L4(T) = 0, L5(T) = 0, P(0) = P0, Sh(0) = Sh0, Sv(0) = Sv0}, fcns, type = numeric);
fcns := {Ih(t), Iv(t), L1(t), L2(t), L3(t), L4(t), L5(t), P(t), 

  Sh(t), Sv(t)}
Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

Hi, i want to plot an ode plot with respect to a parameter (E) in range [0,4] instead of time.
I start with unassign variable E and tried to change 't' for 'E' directly, but its give me error 
did I use the incorrect syntax/ step?
 

restart; with(linalg); with(VectorCalculus); with(DEtools); with(plots)

r := .9; K := 10; beta := .5; g := 0.3e-1; alpha[1] := .4; alpha[2] := 2.2; alpha[3] := 4; d[1] := .1; d[2] := .1; q := 1; E := 2.5; s := .46; delta := 1.5; b := 1.2; p := 1; c := 0.1e-1; mu := 0.25e-1; P0 := 4; M0 := 5; N0 := 3; L0 := 2; T := 20

> dP := VectorCalculus[`+`](VectorCalculus[`+`](VectorCalculus[`*`](VectorCalculus[`*`](r, P(t)), VectorCalculus[`+`](1, VectorCalculus[`-`](VectorCalculus[`*`](P(t), 1/K)))), VectorCalculus[`*`](g, P(t))), VectorCalculus[`-`](VectorCalculus[`*`](beta, P(t)))); dM := VectorCalculus[`+`](VectorCalculus[`+`](VectorCalculus[`*`](beta, P(t)), VectorCalculus[`-`](VectorCalculus[`*`](VectorCalculus[`*`](q, E), M(t)))), VectorCalculus[`-`](VectorCalculus[`*`](d[1], M(t)))); dN := VectorCalculus[`+`](VectorCalculus[`-`](VectorCalculus[`*`](s, N(t))), VectorCalculus[`*`](VectorCalculus[`*`](VectorCalculus[`*`](delta, N(t)), M(t)), 1/VectorCalculus[`+`](M(t), VectorCalculus[`*`](b, N(t))))); dL := VectorCalculus[`+`](VectorCalculus[`*`](VectorCalculus[`+`](alpha[1], VectorCalculus[`-`](VectorCalculus[`*`](VectorCalculus[`*`](alpha[2], L(t)), 1/VectorCalculus[`+`](alpha[3], M(t))))), L(t)), VectorCalculus[`-`](VectorCalculus[`*`](d[2], L(t))));

> satu := diff(P(t), t) = dP; dua := diff(M(t), t) = dM; tiga := diff(N(t), t) = dN; empat := diff(L(t), t) = dL;

> pdb := satu, dua, tiga, empat; fcns := {L(t), M(t), N(t), P(t)};
> Q := dsolve({pdb, L(0) = L0, M(0) = M0, N(0) = N0, P(0) = P0}, fcns, type = numeric, method = rkf45, maxfun = 500000);
> odeplot(Q, [[t, P(t), color = blue], [t, M(t), color = green], [t, N(t), color = red], [t, L(t), color = gold]], t = 0 .. T, numpoints = 100000, thickness = 2); 

its work fine if i plot with respect to time (t), but when i tried to change it for a parameter like E its doesnt work

unassign('E'); 
odeplot(Q, [[E, P(t), color = blue]], E = 0 .. 4, numpoints = 100000, thickness = 2)

Error, (in plots/odeplot) curve is not fully specified in terms of the ODE solution, found additional unknowns {E}
following the output that I expected

attached: captive_breeding.mw

im trying to input this variable with value structure

> xxx := x[2, 0], x[2, 1], x[2, 2], x[1, 0], x[1, 1]

check the type

> whattype(xxx);
exprseq

but when i try to use Get tools from maplets package, its give me error

> xxx:=Get`('xxx1'::exprseq)
invalid argument(s): xxx1::exprseq

is it possible to passing exprseq from maplet input into procedure?
or is there another way to input it (xxx) ?

You are please to check maple file:    robust_rev1.mw

case close, many thanks @tomleslie and @vv

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. Everything works normally and I get a numerical solution for the system

restart;
with(linalg);
with(DEtools);
with(plots);

Model declaration

dx[1] := -T*V*k*u+B*T*r-T*m1+A; dx[2] := T*V*k*u-Ti*m2; dx[3] := N*Ti*m2-V*m3; H := A*T-(1-u)^2+add(dx[i]*L[i], i = 1 .. 3); satu := -(diff(H, T)); dua := -(diff(H, Ti)); tiga := -(diff(H, V)); empat := diff(H, L[1]); lima := diff(H, L[2]); enam := diff(H, L[3])

eq1 := diff(L1(t), t) = -A-(-V(t)*k*u(t)+B*r-m1)*L1(t)-u(t)*k*V(t)*L2(t); eq2 := diff(L2(t), t) = -N*m2*L3(t)+m2*L2(t); eq3 := diff(L3(t), t) = T(t)*k*u(t)*L1(t)-T(t)*k*u(t)*L2(t)+m3*L3(t); eq4 := diff(T(t), t) = -T(t)*V(t)*k*u(t)+B*T(t)*r-T(t)*m1+A; eq5 := diff(Ti(t), t) = T(t)*V(t)*k*u(t)-Ti(t)*m2; eq6 := diff(V(t), t) = N*Ti(t)*m2-V(t)*m3

Value for parameter

 u(t):=-1/(2)*L1(t)*k*V(t)*T(t)+1/(2)*L2(t)*k*V(t)*T(t)+1; s:=10;  m1:=0.02;  m2:=0.5;  m3:=4.4;  r:=0.03;  Tm:=1500;  k:=0.000024;  N:=300;    A:=1;  B:=1-(T(t)+Ti(t))/(Tm): 

Numerical Solution with BVP
 

fcns := {L1(t), L2(t), L3(t), T(t), Ti(t), V(t)}; a := dsolve({eq1, eq2, eq3, eq4, eq5, eq6, L1(20) = 0, L2(20) = 0, L3(20) = 0, T(0) = 800, Ti(0) = 0.4e-1, V(0) = 1.5}, fcns, type = numeric, method = bvp[midrich])

The plot

odeplot(a, [[t, u(t)], [t, V(t)]], 0 .. 20, numpoints = 1000)

Output

Blue: V(t)
Red: u(t)
Based on graphs v (t) and u (t) have negative values, whereas in fact v (t) shows many viruses where it is never negative (this is irrational) and u (t) is limited in interval [0,1]. My question are:
How to provide positive assumptions for the system solution?
So v (t), T (t), Ti (t) are never negative?