restart:

#
# Define gamma as local (don't like doing this!)
#
  local gamma:
#
# Replaced 'indexed' parameters with 'inert subscript'
# parameters - otherwise one gets a problem defining
# both the unindexed 'phi' and the indexed phi[c]
#
  b := .45: mu__t := 1.112: delta := .181: rho := 0.2e-1:
  beta__k := .123: sigma := 0.9e-1: alpha := 0.312e-1: gamma := 0.14e-2:
  phi := .24: xi := .134: A__h := .576: k__1 := 0.1e-2:
  beta__c := 0.1e-1: mu__c := 0.19e-2:   eta := .557: z := .636:
  phi__c := 0.57e-1:
#
# D() is Maple's differential operator replated D(T)
# with DD(T) in the following to avoid confusion
#
  ODE1 := diff(B(T), T) = rho*b-mu__t*B(T)-delta*B(T):
  ODE2 := diff(C(T), T) = delta*B(T)-mu__t*C(T)+1-rho*b-beta__k*C(T)*H(T)+sigma*C(T):
  ODE3 := diff(DD(T), T) = beta__k*C(T)*H(T)-(mu__t+alpha+gamma)*C(T)-phi*C(T):
  ODE4 := diff(E(T), T) = alpha*DD(T)-xi*E(T)-mu__t*E(T):
  ODE5 := diff(F(T), T) = phi*DD(T)+xi*E(T)-mu__t*F(T)-sigma*C(T):
  ODE6 := diff(G(T), T) = (1-k__1)*A__h-beta__c*G(T)*H(T)-mu__c*H(T)-eta*J(T):
  ODE7 := diff(H(T), T) = k__1*A__h+beta__c*G(T)*H(T)-mu__c*H(T)-z*H(T):
  ODE8 := diff(J(T), T) = z*H(T)-(phi__c+eta)*J(T):

  B0 := 100: C0 := 90: D0 := 45: E0 := 38:
  F0 := 100: G0 := 45: H0 := 50: J0 := 20:
#
# Solve system
#
  ans := dsolve( { ODE1, ODE2, ODE3, ODE4, ODE5, ODE6, ODE7, ODE8,
                   B(0) = B0, C(0) = C0, DD(0) = D0, E(0) = E0,
                   F(0) = F0, G(0) = G0, H(0) = H0, J(0) = J0
                 },
                 numeric
               );
#
# Plot solutions for a few of the dependent variablss
# just to show everything is working (more-or-less!)
#
  plots:-odeplot( ans, [T, DD(T)], T=0..5);
  plots:-odeplot( ans, [T, H(T)], T=0..5);
  plots:-odeplot( ans, [T, J(T)], T=0..5);

  restart;
  with(plots):
  eqs1:= diff(r[1](t), t, t)+.68041478*(diff(r[1](t), t))+732.61274*r[1](t)-0.11751824e-2*V(t) = -0.52084940e-2*q(t),
         2.03*10^(-8)*(diff(V(t), t))+4.0650406*10^(-7)*V(t)+0.11751824e-2*(diff(r[1](t), t)) = 0,
         diff(q(t), t, t)+(9.5199780*(q(t)^2-1))*(diff(q(t), t))+1006.9998*q(t) = -2192.5091*(diff(r[1](t), t, t));

  ics1:= V(0) = 0, q(0) = 0.1e-2, (D(q))(0) = 0, r[1](0) = 0, (D(r[1]))(0) = 0;\
  sol1:= dsolve( [eqs1, ics1], numeric);
#
# Plot results separately, just because r[1](t)
# is very much smaller than V(t) and q(t). In
# a "combined" plot it "looks" as if r[1](t)~0
#
  odeplot( sol1,
           [ t, r[1](t) ],
           t=0..2,
           axes=boxed,
           title="r1(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );
  odeplot( sol1,
           [ t, V(t) ],
           t=0..2,
           axes=boxed,
           title="V(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );
  odeplot( sol1,
           [ t, q(t) ],
           t=0..2,
           axes=boxed,
           title="q(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );

#
# This group fails because there are too few equations
# (three) for the number of dependent variables
# (four), as well as the fact that there are no
# initial conditions for the dependent variable
# r[1](t)
#    
  eqs2:=diff(r[2](t), t, t)+3.3737665*(diff(r[2](t), t))+52672.540*r[2](t)-0.93177795e-2*V(t) = 0.18358493e-2*q(t),
        2.03*10^(-8)*(diff(V(t), t))+4.0650406*10^(-7)*V(t)+0.11751824e-2*(diff(r[1](t), t))+0.93177795e-2*(diff(r[2](t), t)) = 0,
        diff(q(t), t, t)+(9.5199780*(q(t)^2-1))*(diff(q(t), t))+1006.9998*q(t) = -2192.5091*(diff(r[1](t), t, t))+772.79848*(diff(r[2](t), t, t));
  ics2 := V(0) = 0, q(0) = 0.1e-2, (D(q))(0) = 0, r[2](0) = 0, (D(r[2]))(0) = 0;
  sol2:=dsolve( [eqs2, ics2], numeric);

Error, (in DEtools/convertsys) ODE system is insufficient to determine values for all dependent variables in the system

#
# Combine the two equation systems to produce six
# equations in four dependent variables, which might(?)
# be solvable if the resulting equation system is
# "consistent"
#
# However Maple cannot convert to an explicit first
# order system, so cannot solve.
#
# This fails - query inconsistency?
#
  sol2:=dsolve( {eqs1, eqs2, ics1, ics2}, numeric);

Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

#
# With no real justification, just select the first
# equation from 'eqs1'. Combine this with 'eqs2' and
# combine the ics
#
# At least this has the "correct" number of equations
# for the number of dependent variables (and the
# correct number of ics)
#
# This works, but is it meaningful
#
  sol2:=dsolve( {eqs1[1], eqs2, ics1, ics2}, numeric);
#
# Plot separately just because r[](t) and r[2](t)
# are very small compared with V(t) and q(t)
#
  odeplot( sol2,
           [ t, r[1](t) ],
           t=0..2,
           axes=boxed,
           title="r1(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );
  odeplot( sol2,
           [ t, r[2](t) ],
           t=0..2,
           axes=boxed,
           title="r2(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );
  odeplot( sol2,
           [ t, V(t) ],
           t=0..2,
           axes=boxed,
           title="V(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );
  odeplot( sol2,
           [ t, q(t) ],
           t=0..2,
           axes=boxed,
           title="q(t) vs t",
           titlefont=[times, bold, 20],
           numpoints=1000
         );

  restart;
  with(Statistics):

  X := RandomVariable(Normal(0, 1));
#
# Probability that 0 < X < 0.5
#
  Probability( { X > 0, X < 0.5 } );
#
# Probability that -2.24 < X < 1.12
#
  Probability( { X > -2.24, X < 1.12 });
#
# With X a random variable, normally distributed,
# unit standard deviation: compute value of x
# such that
#
#    prob(-x<X<x)=0.847
#
  f:=x->Probability( { X > -x, X < x } )-0.847:
  fsolve(f);

#
# Doing everything symbolically. NB if my interpretation
# of BCs, ICs is correct, then solution only exists in the
# range x=0..2
#
  restart;
  pde:=diff(u(x,t),t,t)=9*diff(u(x,t),x,x):
  bcs:=u(0,t)=0, u(2,t)=0:
  ics:= u(x,0) = piecewise( x<0,  undefined,
                            x<1,  x*(2-x),
                            x<=2, x,
                            undefined
                          ),
        D[2](u)(x,0) = 0:
  sol:=pdsolve( [pde, bcs, ics]);

#
# Since the soluton is an infinite sum with no
# closed form, this solution is "difficult" to
# plot. However one can easily extract the summand,
# and easily plot the first few terms of the sum
# ie n=1, n=2, n=3, n=4....
#
  solfn:=unapply( op([2,1],sol), [n,x,t]);
  nterms:= 4:
  for j from 1 by 1 to nterms do
      plot3d( solfn(j, x, t),
              x = 0..2,
              t = 0..10,
              title = typeset("n=", j," term"),
              titlefont = [times, bold, 20]
            )
  od;
#
# Alternatively, one could produce a plot of the sum
# of the first few terms
#
  nterms:= 50;
  plot3d( add
          ( solfn(j,x,t),
            j=1..nterms
          ),
          x = 0..2,
          t = 0..10,
          title = typeset("Sum of first ", nterms," terms"),
          titlefont = [times, bold, 20]
        );
#
# Or an animation, in time
#
  plots:-animate( plot,
           [ add
             ( solfn(j,x,t),
               j=1..nterms
             ),
             x = 0..2
           ],
           t = 0..10,
           frames=100
        );

##########################################
# Doing everything numerically
#
  restart;
  np:=1000:
  pde:=diff(u(x,t),t,t)=9*diff(u(x,t),x,x):
  bcs:=u(0,t)=0, u(2,t)=0:
  ics:= u(x,0) = piecewise( x<1,  x*(2-x),
                            x<=2, x
                          ),
        D[2](u)(x,0) = 0:
  
  sol:= pdsolve( pde, [bcs, ics], numeric, time=t, range=0..2, timestep=0.01, spacestep=5/np);
  sol:-plot3d( t=0..10, x=0..2, numpoints=10000);
#
# Alternatively one can produce an animation in time
# using something like
#
  sol:-animate( t=10,
                frames=100,
                title="time = %f"
              );

#
# Doing everything symbolically. NB if my interpretation
# of BCs, ICs is correct, then solution only exists in the
# range x=0..5
#
  restart;
  pde:=diff(u(x,t),t,t)=4*diff(u(x,t),x,x):
  bcs:=u(0,t)=0, u(5,t)=0:
  ics:= u(x,0) = piecewise( x<0, undefined,
                            x<4, 0,
                            x<=5,  5-x,
                            undefined
                          ),
        D[2](u)(x,0) = 0:
  sol:=pdsolve( [pde, bcs, ics]);

#
# Since the soluton is an infinite sum with no
# closed form, this solution is "difficult" to
# plot. However one can easily extract the summand,
# and easily plot the first few terms of the sum
# ie n=1, n=2, n=3, n=4....
#
  solfn:=unapply( op([2,1],sol), [n,x,t]);
  nterms:= 4:
  for j from 1 by 1 to nterms do
      plot3d( solfn(j, x, t),
              x = 0..5,
              t = 0..10,
              title = typeset("n=", j," term"),
              titlefont = [times, bold, 20]
            )
  od;
#
# Alternatively, one could produce a plot of the sum
# of the first few terms
#
  nterms:= 50;
  plot3d( add
          ( solfn(j,x,t),
            j=1..nterms
          ),
          x = 0..5,
          t = 0..10,
          title = typeset("Sum of first ", nterms," terms"),
          titlefont = [times, bold, 20]
        );
#
# Or an animation, in time
#
  plots:-animate( plot,
           [ add
             ( solfn(j,x,t),
               j=1..nterms
             ),
             x = 0..5
           ],
           t = 0..10,
           frames=100
        );

##########################################
# Doing everything numerically
#
  restart;
  np:=1000:
  pde:= diff(u(x,t),t,t)=4*diff(u(x,t),x,x):
  bcs:= u(0,t)=0, u(5,t)=0:
  ics:= u(x,0) = piecewise( x<4, 0,
                            x<=5,  5-x
                          ),
        D[2](u)(x,0) = 0:
  sol:= pdsolve( pde, [bcs, ics], numeric, time=t, range=0..5, timestep=0.01, spacestep=5/np);
  sol:-plot3d( t=0..10, x=0..5, numpoints=10000);
#
# Alternatively one can produce an animation in time
# using something like
#
  sol:-animate( t=10,
                frames=100,
                title="time = %f"
              );

  restart;
  with(VectorCalculus):
  vs1 := VectorSpace('cartesian', [0,0]):
  vs2 := VectorSpace('cartesian', [1,2]):
#
# To make all the vectors the "same" width
# use the 'width' option, which will be
# applied to all vectors in the list
#
  PlotVector( [ vs1:-Vector([1,2]),
                vs1:-Vector([3,2]),
                vs2:-Vector([2,0])
              ],
              color=[red, green, blue],
              width=[0.5, relative=false]
            );
#
# If you want to specify widths individually
# then you have to use multiple PlotVector()
# commands, with  individual 'width'
# specifications
#
  p1:= PlotVector( vs1:-Vector([1,2]),
                   color=red,
                   width=[0.5, relative=false]
                 ):
  p2:= PlotVector( vs1:-Vector([3,2]),
                   color=green,
                   width=[0.25, relative=false]
                 ):
  p3:= PlotVector( vs2:-Vector([2,0]),
                   color=blue,
                   width=[0.1, relative=false]
                 ):
  plots:-display( [p1, p2, p3]);

  restart;

#
# Define parameters: set up and solve the ODE
# system
#
  a:= .5; b:= 0.3e-1; c:= .5; d:= 0.1e-1:
  ode1:= diff(x(t), t) = a*x(t)-b*x(t)*y(t):
  ode2:= diff(y(t), t) = -c*y(t)+d*x(t)*y(t):
  sol:= dsolve([ode1, ode2, x(0) = 50, y(0) = 10], numeric, output=listprocedure);
#
# Notice the output format of the above. One
# can now assign the individual solution
# procedures to names for ease of access,
# as in the following
#
  xx:= eval( x(t), sol);
  yy:= eval( y(t), sol);

#
# Now we can use these simple procedures in the
# same way any other procedures might be used.
#
# So if you just want the value of x(t) or y(t)
# at any supplied times
#
  xx(20);
  xx(25);
  yy(12);
  yy(33);

#
# Or, for example, plotting without using
# ODEplot
#
  plot( [xx(t), yy(t)],
        t=0..50,
        color=[red, blue]
      );
#
# Of course you can still use odeplot if you want,
# but now it is looking like the *harder* way to
# do it!
#
  plots:-odeplot( sol,
                  [ [t, x(t)],
                    [t, y(t)]
                  ],
                  t=0..50,
                  color=[red, blue]
                )

#
# Or finding specific values within a given range
#
  fsolve( xx(t)=60, t=0..10);
  fsolve( yy(t)=20, t=0..10);

#
# Or finding the first 10 time-values when x(t)=60
#
  with(RootFinding):
  tstart:=0:
  for i from 1 to 10 do
      tstart:=NextZero( t->xx(t)-60, tstart);
  od;

restart;
pde := diff(u(r, theta, t), t) = (diff(u(r, theta, t), r)+r*(diff(u(r, theta, t), r, r))+(diff(u(r, theta, t), theta, theta))/r)/(25*r);
bc_on_r := D[1]*u(1,theta,t)=0;
bc_on_theta:= u(r,0,t)=0, u(r,Pi,t)=0;
ic := u(r,theta,0)=(r-1/3*r^3)*sin(theta);
sol:=pdsolve([pde, bc_on_r,bc_on_theta,ic], u(r, theta, t), HINT = boundedseries(r = [0]))

restart;
pde := diff(u(r, theta, t), t) = (diff(u(r, theta, t), r)+r*(diff(u(r, theta, t), r, r))+(diff(u(r, theta, t), theta, theta))/r)/(25*r);
bc_on_r := u(1,theta,t)=0;
bc_on_theta:= u(r,0,t)=0, u(r,Pi,t)=0;
ic := u(r,theta,0)=(r-1/3*r^3)*sin(theta);
sol:=pdsolve([pde, bc_on_r,bc_on_theta,ic], u(r, theta, t), HINT = boundedseries(r = [0]))

restart;
pde := diff(u(r, theta, t), t) = (diff(u(r, theta, t), r)+r*(diff(u(r, theta, t), r, r))+(diff(u(r, theta, t), theta, theta))/r)/(25*r);
bc_on_r := D[1](u)(1,theta,t)=0;
bc_on_theta:= u(r,0,t)=0, u(r,Pi,t)=0;
ic := u(r,theta,0)=(r-1/3*r^3)*sin(theta);
sol:=pdsolve([pde, bc_on_r,bc_on_theta,ic], u(r, theta, t), HINT = boundedseries(r = [0]))

#
# Initialize and set up necessary functions
#
  restart;
  with(Optimization):
  y := (p, e) -> alpha-beta*p+k*e:
  g := e-> (1/2)*mu*e^2:
  Lambda:= z-> int((z-u)*phi(u), u=0..z):
  PI := unapply( (p-c)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p,e,z]):
  Pi_rd:= unapply( (p-w)*(y(p, e)-z)-(p-v)*Lambda(z)-g(e), [p, e, z]):

#
# Set up k-values and optimizer conditions
#
  kVals:=[1, 2, 3]:
  optConds:=  p = 0 .. 60,
              e = -1000 .. 30,
              z = -10000 .. 25,
              initialpoint = {e = 0, p = 10, z = 0},
              maximize:

  for ind from 1 by 1 to numelems(kVals) do
    #
    # Set up parameters
    #
      pars:= [ k = kVals[ind],
               beta = 2,
               alpha = 50,
               mu = 10,
               c = 5,
               v = 1,
               phi(u) = 1/2,
               w=10
             ]:
      Int_Profit[ind]:= NLPSolve
                        ( eval
                          ( PI(p, e, z),
                            pars
                          ),
                          optConds
                       );
      qopt[ind]:= eval(z+y(p, e), [ Int_Profit[ind][2][], pars[]]);
      Re_Profit[ind]:= NLPSolve
                       ( eval
                         ( Pi_rd(p, e, z),
                           pars
                         ),
                         optConds
                       );
      qd_opt[ind]:= eval(z+y(p, e), [ Re_Profit[ind][2][], pars[]]);
      M_Profit[ind]:= eval((w-c)*qd_opt[ind], pars);
      Whole_Profit[ind]:= M_Profit[ind]+Re_Profit[ind][1];
      t[ind]:= eval((Int_Profit[ind][1]-Re_Profit[ind][1])/(w-c), pars);
  od:

#
# How to access results, eg for Re_Profit.
#
# Others ie Int_Profit
#           Whole_Profit
#           qopt
#           qd_opt
#           t
#
# are accessed similarly
#
  Re_Profit[1];
  Re_Profit[2];
  Re_Profit[3];

#
# Set some (more-or-less-random) values
# for parameters
#
 K:=1: k:=1: lambda:=1: mu:=1: kappa:=1:

plot3d(kappa*sech(kappa*(-k^2*t+x)), x = -10 .. 10, t = -10 .. 10);
plot3d((2*(-K^2+k^2))*(K*cosh(k*(-k^2*t+x))+k*cosh(K*(-K^2*t+x)))/((K+k)^2*cosh(k*(-k^2*t+x)-K*(-K^2*t+x))+(-K+k)^2*cosh(k*(-k^2*t+x)+K*(-K^2*t+x))+4*K*k), x = -10 .. 10, t = -10 .. 10);
plot3d(2*(diff(arctan(lambda*sin(mu*((-3*lambda^2+mu^2)*t+x))*sech(lambda*((-lambda^2+3*mu^2)*t+x))/mu), x)), x = -10 .. 10, t = -10 .. 10);

  restart;
  ode:=diff(y(t),t)+3*y(t)=1;
  ic:= y(0)=0;
#
# Solve the problem numerically
#
  soln:=dsolve( [ode, ic], numeric, output=listprocedure):
  yt:=eval( y(t), soln):
#
# Solve the problem analytically
#
  sola:= dsolve([ode, ic]):
#
# Plot the analytic and numeric solution on
# the same graph - can't tell the difference
# and include the OP's (completely wrong)
# "solution" whihc is colored green in the
# following
#
  plot( [yt(t), rhs(sola), -3*t^2/2+t],
         t=0..1,
         color=[red, blue, green]);
#
# So plot the difference between the analytic
# numeric solution on the same graph (and be
# sure to check the vertical scale before
# complaining!)
#
  plot( yt(t)-rhs(sola), t=0..2);