Here I attached the Code, Thank you for your assistance.

Download Simulation_Code.mw

Dear @Kitonum,

May you halp me solving may error. I itied however, I cant see the problem.

Thank you in advance ,

##############################################

restart;
with(plots);
alpha[1] := .4; alpha[2] := .4; rho[2] := .4; PI := 50; delta[2] := .27; `ϑ`[2] := .27; sigma[1] := .15; sigma[2] := .285; sigma[3] := 0.35e-1; phi[vl] := .18; Lambda := 50;
beta[1] := .15; rho[1] := .2; delta[1] := .15; `ϑ`[1] := .1468; `ϕ`[vd] := .71; Omega := 50; mu := 0.128e-1; tau := 0.25e-1; epsilon := .1; eta := 0.2e-1; k := 0.2e-2;
a[1] := 1; a[2] := 1; a[3] := 2; a[4] := 2; a[5] := 2; w[1] := 100; w[2] := 100; w[3] := 100; w[4] := 200; w[5] := 300; T := 3;

u[1] := min(max(0, q), 1)*cont; q := (c*beta[1]*i[vd](t)*(s(t)*(lambda[2](t)-lambda[1](t))+i[l](t)*(lambda[4](t)-lambda[3](t)))+(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))*((lambda[3](t)-lambda[1](t))*s(t)+(lambda[4](t)-lambda[2](t))*i[d](t))+c*rho[1]*i[d](t)*s[vd](t)*(lambda[9](t)-lambda[8](t))+rho[2]*i[l](t)*(s[vl](t))(lambda[11](t)-lambda[9](t)))/w[1]; u[2] := min(max(0, v), 1); v := (i[d](t)*(lambda[2](t)-lambda[5](t))+(i[dl](t))(lambda[4](t)-lambda[5](t)-lambda[7](t)))/w[2]; u[3] := min(max(0, z), 1); z*`≔`*(i[l](t)*(lambda[3](t)-lambda[6](t))+(i[dl](t))(lambda[4](t)-lambda[6](t)-lambda[7](t)))/w[3]; u[4] := min(max(0, n), 1); n := (lambda[8](t)*s[vd](t)+lambda[9](t)*i[vd](t))/w[4]; u[5] := min(max(0, p), 1); p := (lambda[10](t)*s[vd](t)+lambda[9](t)*i[vd](t))/w[5];

u[1] := 0; u[2] := 0; u[3] := 0; u[5] := 0; u[4] := 0;
NULL;
sys := diff(s(t), t) = PI+sigma[1]*r[d](t)+sigma[2]*r[l](t)+sigma[3]*r[dl](t)-(1-u[1])*c*beta[1]*i[vd](t)*s(t)-(1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))*s(t)-mu*s(t), diff(i[d](t), t) = (1-u[1])*c*beta[1]*i[vd](t)*s(t)-(1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))*i[d](t)-(`ϑ`[1]+u[3])*i[d](t)-(delta[1]+mu)*i[d](t), diff(i[l](t), t) = (1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))*s(t)-(1-u[1])*c*beta[1]*i[vd](t)*i[l](t)-(`ϑ`[2]+u[3])*i[l](t)-(delta[2]+mu)*i[l](t), diff(i[dl](t), t) = (1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))*i[d](t)+(1-u[1])*c*beta[1]*i[vd](t)*i[l](t)-(u[2]+u[3]+epsilon)*i[dl](t)-(delta[1]+delta[2]+mu)*i[dl](t), diff(r[d](t), t) = (`ϑ`[1]+u[2])*i[d](t)+(epsilon*eta+u[2])*i[dl](t)-(sigma[1]+mu)*r[d](t), diff(r[l](t), t) = (`ϑ`[2]+u[3])*i[l](t)+(epsilon*k*(1-eta)+u[3])*i[dl](t)-(sigma[2]+mu)*r[l](t), diff(r[dl](t), t) = (epsilon*(k-1)*(1-eta)+u[2]+u[3])*i[dl](t)-(sigma[3]+mu)*r[dl](t), diff(s[vd](t), t) = Omega-(1-u[1])*tau*rho[1]*i[d](t)*s[vd](t)-(u[4]+`ϕ`[vd])*s[vd](t), diff(i[vd](t), t) = (1-u[1])*tau*rho[1]*i[d](t)*s[vd](t)-(u[4]+`ϕ`[vd](t))*i[vd](t), diff(s[vl](t), t) = Lambda-(1-u[1])*rho[2]*i[l](t)*s[vl](t)-(u[5]+phi[vl])*s[vl](t), diff(i[vl](t), t) = (1-u[1])*rho[2]*i[l](t)*s[vl](t)-(u[5]+phi[vl])*i[vl](t), diff(lambda[1](t), t) = lambda[1](t)*((1-u[1])*tau*beta[1]*i[vd](t)+(1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))+mu)-lambda[2](t)*(1-u[1])*c*beta[1]*i[vd](t)-lambda[3](t)*(1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t)), diff(lambda[2](t), t) = -a[1]+lambda[2](t)*((1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))+`ϑ`[1]+u[2]+delta[1]+mu)-lambda[4](t)*(1-u[1])*(alpha[1]*(i[l](t)+i[dl](t))+alpha[2]*i[vl](t))-lambda[5](t)*(`ϑ`[1]+u[2])+lambda[8](t)*(1-u[1])*c*rho[1]*s[vd](t)-lambda[9](t)*(1-u[1])*c*rho[1]*s[vd](t), diff(lambda[3](t), t) = -a[2]+lambda[1](t)*(1-u[1])*alpha[1]*s(t)+lambda[2](t)*(1-u[1])*alpha[1]*i[d](t)-lambda[3](t)*((1-u[1])*alpha[1]*s(t)-(1-u[1])*c*beta[1]*i[vd](t)-`ϑ`[2]-u[3]-delta[2]-mu)-lambda[4](t)*((1-u[1])*alpha[1]*i[d](t)+(1-u[1])*c*beta[1]*i[vd](t))-lambda[6](t)*(`ϑ`[2]+u[3])+lambda[10](t)*(1-u[1])*rho[2]*s[vl](t)-lambda[11](t)*(1-u[1])*rho[2]*s[vl](t), diff(lambda[4](t), t) = -a[3]+lambda[1](t)*(1-u[1])*alpha[1]*s(t)+lambda[2](t)*(1-u[1])*alpha[1]*i[d](t)-lambda[3](t)*(1-u[1])*alpha[1]*s(t)-lambda[4](t)*((1-u[1])*alpha[1]*i[d](t)-u[2]-u[3]-epsilon-delta[1]-delta[2]-mu)-lambda[5](t)*(epsilon*eta+u[2])-lambda[6](t)*(epsilon*k*(1-eta)+u[3])-lambda[7](t)*(epsilon*(k-1)*(1-eta)+u[2]+u[3]), diff(lambda[5](t), t) = -lambda[1](t)*sigma[1]+lambda[5](t)*(sigma[1]+mu), diff(lambda[6](t), t) = -lambda[1](t)*sigma[2]+lambda[6](t)*(sigma[2]+mu), diff(lambda[7](t), t) = -lambda[1](t)*sigma[3]+lambda[7](t)*(sigma[3]+mu), diff(lambda[8](t), t) = lambda[8](t)*((1-u[1])*c*rho[1]*i[d](t)+u[4]+`ϕ`[vd])-lambda[9](t)*(1-u[1])*c*rho[1]*i[d](t), diff(lambda[9](t), t) = -a[4]+(lambda[1](t)-lambda[2](t))*(1-u[1])*tau*beta[1]*s(t)+(lambda[3](t)-lambda[4](t))*(1-u[1])*tau*beta[1]*i[l](t)+lambda[9](t)*(u[4]+`ϕ`[vd]), diff(lambda[10](t), t) = lambda[10](t)*((1-u[1])*rho[2]*i[l](t)+u[5]+phi[vl])-lambda[11](t)*(1-u[1])*rho[2]*i[l](t), diff(lambda[11](t), t) = -a[5]+(lambda[1](t)-lambda[3](t))*(1-u[1])*alpha[2]*s(t)+(lambda[2](t)-lambda[4](t))*(1-u[2])*alpha[2]*i[d](t)+lambda[11](t)*(u[5]+phi[vl]), s(0) = 10000, i[d](0) = 10, r[d](0) = 10, s[vd](0) = 500, i[vd](0) = 100, i[dl](0) = 100, i[l](0) = 100, r[l](0) = 10, r[dl](0) = 10, s[vl](0) = 500, i[vl](0) = 100, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0, lambda[4](T) = 0, lambda[5](T) = 0, lambda[6](T) = 0, lambda[7](T) = 0, lambda[8](T) = 0, lambda[9](T) = 0, lambda[10](T) = 0, lambda[11](T) = 0:

p1 := dsolve({sys}, type = numeric, method = bvp[midrich], abserr = .4, continuation = cont, maxmesh = 2400):
p2o := odeplot(p1, [t, e(t)], 0 .. 4, numpoints = 100, labels = ["Time (months)", "Carrier Human"], labeldirections = [horizontal, vertical], style = line, color = red, axes = boxed);

     tomleslie 5927

Thank you for the brief, but I have a question again

You mean it cannot be solved if a recursive assignment exists?

   @ tomleslie 5927

Thank you too mach but I cant understand what you want to show me. May you make it clear.

If possible send me after putting the your comments. Thank you again

Dear acer 17804  and tomleslie 5927

I would like to thank for the comment you gave me. However, I did all the corrections but it is still giving the same error.

May I got your help please?

restart;
with(plots);
beta := 0.143e-1; delta := 0.7e-1; PI := 2.5; mu := 0.16e-1; tau := .7; rho := 0.4e-3; epsilon := .15; sigma[2] := 0.2e-1; sigma[1] := 0.1e-3; eta := 0.1e-1;
a[1] := 2; a[2] := 5; w[1] := 10; w[2] := 50; w[3] := 20; T := 4;

u[1] := min(max(0, z), 1); z := beta*(1-u[3])*s(t)*(sigma[2]*e(t)+sigma[1]*i(t))*(lambda[2](t)-lambda[1](t))/w[1]; u[2] := min(max(0, c), 1); c := (lambda[2](t)*e(t)+lambda[4](t)*((1-tau)*delta*e(t)+epsilon*i(t))+lambda[3](t)*(delta*tau*e(t)+epsilon*i(t)))/w[2]; u[3] := min(max(0, p), 1); p := beta*(1-u[1])*s(t)*(sigma[2]*e(t)+sigma[1]*i(t))*(lambda[2](t)-lambda[1](t))/w[3];
Error, recursive assignment
NULL;
sys := diff(s(t), t) = PI+eta*r(t)-((1-u[1])*(1-u[1]))*beta*(sigma[1]*i(t)+sigma[2]*e(t))*s(t)-mu*s(t), diff(e(t), t) = ((1-u[1])*(1-u[1]))*beta*(sigma[1]*i(t)+sigma[2]*e(t))*s(t)-(u[2]+u[3]+delta+mu)*e(t), diff(i(t), t) = (1-u[2])*tau*delta*e(t)-(u[2]+u[3]+epsilon+rho+mu)*i(t), diff(r(t), t) = (1-u[2])*(1-tau)*delta*e(t)+(1-u[2])*epsilon*i(t)-(mu+eta)*r(t), diff(lambda[1](t), t) = lambda[1](t)*((1-u[3])*(1-u[1])*beta*(sigma[1]*i(t)+sigma[2]*e(t))-mu)-lambda[2](t)*(1-u[3])*(1-u[1])*beta*(sigma[1]*i(t)+sigma[2]*e(t)), diff(lambda[2](t), t) = -a[1]+lambda[1](t)*(1-u[3])*(1-u[1])*beta*sigma[2]*s(t)-lambda[2](t)*((1-u[3])*(1-u[1])*beta*sigma[2]*s(t)-u[2]-u[3]-delta-mu)-lambda[3](t)*(1-u[2])*tau*delta-lambda[4](t)*(1-u[2])*(1-tau)*delta, diff(lambda[3](t), t) = -a[2]+lambda[1](t)*(1-u[3])*(1-u[1])*beta*sigma[1]*s(t)-lambda[2](t)*(1-u[3])*(1-u[1])*beta*sigma[1]*s(t)+lambda[3](t)*(u[2]+u[3]+epsilon+rho+mu)-lambda[4](t)*(1-u[2])*epsilon, diff(lambda[4](t), t) = -lambda[1](t)*eta+lambda[4](t)*(mu+eta), s(0) = 1000000, e(0) = 1000, i(0) = 89, r(0) = 16, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0, lambda[4](T) = 0;
p1 := dsolve({sys}, type = numeric, method = bvp[midrich], abserr = .1);
Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 9, got 8

Error, (in dsolve/numeric/bvp/convertsys) e(t) and e cannot both appear in the given ODE
p2o := odeplot(p1, [t, i(t)], 0 .. 2, numpoints = 100, labels = ["Time (months)", "Infectives"], labeldirections = [horizontal, vertical], style = line, color = red, axes = boxed);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
plots[display](p2o);

Dear dharr 1043 ,

Thank you for your inspection, It is working now.

  Dear  Rouben Rostamian 4135

Thank you for your comment, I mean to plot such like Figures,

@dharr 679,

Thank you for your help. As you said I try all what you have commented however can't get a solution. 

would you please give me other possible comments.

Regards,

@dharr 654,

This is to attach the worksheet Infection_model.mw
 

Download Infection_model.mw
 

Download Infection_model.mw
 

Download Infection_model.mw

Dear  Kitonum 12774 ,

Thank you 

Now it says Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging.