restart:

  interface(rtablesize=10):

#
# Define gamma as local (don't like doing this!)
#
  local gamma:local pi:

#
# Replaced 'indexed' parameters with 'inert subscript'
# parameters - otherwise one gets a problem defining
# both the unindexed 'phi' and the indexed phi[c]
#

if false then
  A__h := .18: beta__1 := 0.8354e-1: psi := 0.7258e-1: mu__r := 0.51e-1: sigma := .165: alpha := .65: xi := 0.5e-1: gamma := .131: A__b := .7241: beta__o := 0.65e-1: mu__b = 0.17e-1:
end if:

#
# D() is Maple's differential operator replated D(T)
# with DD(T) in the following to avoid confusion
#
  ODE1 := diff(S(T), T) = A__h-psi*beta__1*S(T)*G(T)-mu__r*S(T):
  ODE2 := diff(G(T), T) = psi*beta__1*S(T)*G(T)-sigma*psi*beta__1*R(T)*B(T)-(alpha+xi+mu__r)*G(T):
  ODE3 := diff(A(T), T) = alpha*G(T)-(gamma+mu__r)*A(T):
  ODE4 := diff(R(T), T) = gamma*A(T)+sigma*psi*beta__1*R(T)*B(T)- mu__r*R(T):
  ODE5 := diff(C(T), T) = A__b - psi*beta__o*C(T)*I(T)- mu__b*C(T):
  ODE6 := diff(B(T), T) = psi*beta__o*C(T)*I(T)- mu__b*B(T):
 

 
if false then
  S0 := 100: G0 := 190: A0 := 45: R0 := 20: C0 := 35: B0 := 25:
end if:

# system + ic


sys := { ODE1, ODE2, ODE3, ODE4, ODE5, ODE6,
                   S(0) = S0, G(0) = G0, A(0) = A0, R(0) = R0, C(0) = C0, B(0) = B0

                 }:

params := convert(indets(sys, name) minus {T}, list);

#
# Solve system
#
  ans := dsolve( { ODE1, ODE2, ODE3, ODE4, ODE5, ODE6,
                   S(0) = S0, G(0) = G0, A(0) = A0, R(0) = R0, C(0) = C0, B(0) = B0
               
                 },
                 parameters = params,
                 numeric
               );

Error, (in dsolve/numeric) 'parameters' must be specified as a list of unique unassigned names

# define parameter values

DefaultValues := [ A__h = .18, beta__1 = 0.8354e-1, psi = 0.7258e-1, mu__r = 0.51e-1, sigma = .165, alpha = .65, xi = 0.5e-1, gamma = .131, A__b = .7241, beta__o = 0.65e-1, mu__b = 0.17e-1,
     S0 = 100, G0 = 190, A0 = 45, R0 = 20, C0 = 35, B0 = 25   
]:

# Form the list of numerical values ordered as params

ValuatedParams := subs(DefaultValues, params)
   

Warning, inserted missing semicolon at end of statement

# Instanciate the solution for this set of values

ans(parameters=ValuatedParams)

Warning, inserted missing semicolon at end of statement

#
# Plot solutions for a few of the dependent variablss
# just to show everything is working (more-or-less!)
#

  plots:-odeplot( ans, [T, S(T)] , T=0..10);
  plots:-odeplot( ans, [T, G(T)] , T=0..10);
  plots:-odeplot( ans, [T, A(T)] , T=0..10);
   plots:-odeplot( ans, [T,R(T)] , T=0..10);
  plots:-odeplot( ans, [T, C(T)] , T=0..10);
  plots:-odeplot( ans, [T, B(T)] , T=0..10);


if false then
  plots:-odeplot( ans, [T, S(T)] , T=0..10);
  plots:-odeplot( ans, [T, G(T)] , T=0..10);
  plots:-odeplot( ans, [T, A(T)], T=0..10);
  plots:-odeplot( ans, [T, R(T)] , T=0..10);
  plots:-odeplot( ans, [T, C(T)] , T=0..10);
  plots:-odeplot( ans, [T, B(T)] , T=0..10);
 


end if:

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

# How two use the parametric solution?
#
# Here an example with only one variable parameter, but this could be generalized to many
# (the main problem is then more a problem of "readability" than a technical one).
#

f := proc(VariableName::symbol, ParamName::symbol, ParamValues::list, params, ParamDefault, sol)
  local V, k, p, MyChoice, ValuatedParams:
  V := unapply(VariableName(T), T):
  k := 0:
  for p in ParamValues do
    k := k+1:
    MyChoice       := map(u -> if lhs(u)=ParamName then lhs(u)=p else u end if, ParamDefault):
    ValuatedParams := subs(MyChoice, params):
    sol(parameters=ValuatedParams):
    plot||k := plots:-odeplot(
                               sol,
                               [T, V(T)] ,
                               T=0..5,
                               legend=p,
                               color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])
                             );
  end do:
  plots:-display(seq(plot||k, k=1..numelems(ParamValues)));
end proc:

f(B, epsilon, [0.0001, 0.334, 5.078], params, DefaultValues, ans);

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

f(B, beta, [0.1, 0.28, 0.35], params, DefaultValues, ans);

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

f(B, alpha , [0.23,0.28, 0.32], params, DefaultValues, ans);

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

f(J, tau , [0.04, 0.05, 0.06], params, DefaultValues, ans);

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

f(C, theta__h , [0.6, 0.7902, 0.8], params, DefaultValues, ans);

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution