M:=Matrix([ lhs~(sol(0)),
            seq(rhs~(sol(j)), j=0..7, 0.01)
          ]
        ):
M[..,3];

  restart:

  with(plots):
  with(plottools):
  with(DEtools):

  eqn1 := diff(V(t), t) = pi*p - (alpha + mu)*V(t),
          V(0) = ic1:

  eqn2 := diff(S(t), t) = alpha*V(t) + (1 - p)*pi - beta*S(t)*In(t)/N - mu*S(t),
          S(0) = ic2:

  eqn3 := diff(In(t), t) = beta*S(t)*In(t)/N - (mu + delta + localgamma)*In(t),
          In(0) = ic3:

  eqn4 := diff(R(t), t) = localgamma*In(t) - mu*R(t),
          R(0) = ic4:

  pi := 487845:
  pVals := [0.2, 0.5, 0.7, 0.8]:
  alpha := 0.054:
  beta := 0.955:
  mu := 0.005:
  delta := 0.03:
  localgamma := 0.935:
  ic1 := 484465:
  ic2 := 31999760:
  ic3 := 26305:
  ic4 := 12470:
  N:=10000:
  cols:=[red, blue, green, black]:
  for k from 1 by 1 to numelems(pVals) do
      p:=pVals[k]:
      dsol := dsolve([eqn1, eqn2, eqn3, eqn4], numeric):
      plt[k]:=odeplot( dsol,
                       [t, rhs(eqn1[1])],
                       t=0..200,
                       color=cols[k]
                     ):
  od:
  display( [seq( plt[j], j=1.. numelems(pVals))],
           legend=[seq( 'p'=pVals[j], j=1..numelems(pVals))],
           labels=[ typeset(t), typeset(lhs(eqn1[1])) ],
           labelfont=[times, bold, 16]
         );

  restart;
  f:= n->expand( (3+2*x+3*x^2)^n );

  coeffs( f(1), x);
  coeffs( f(2), x);

restart:

  Cond:= U(0) = lambda*(D(U))(0), Theta(0) = A+g*(D(Theta))(0),
              U(1) = 0, Theta(1) = B:

OdeSys1:= ((diff(V(Y),Y,Y))/(a6*a7)+Phi(Y)+N*(Phi(Y)*Phi(Y))-a8*(M*M)*V(Y)/a7-(kp*kp)*V(Y)/(a6*a7)),

              (a10*(diff(Phi(Y),Y,Y))/a9+Pr*Ec*((diff(V(Y),Y))*(diff(V(Y),Y))+(V(Y)*V(Y))*(kp*kp))/(a6*a7))
            :

  for j to numelems(MVals) do
      Ans[j] := dsolve( eval([OdeSys, Cond], M = MVals[j]),
                        numeric,
                        output = listprocedure
                      );
      Theta_b[j] := evalf
                    ( Int
                      ( eval
                        ( U(Y)*Theta(Y),
                          Ans[j]
                        )(Y),
                        Y = 0 .. 1
                      )
                      /
                      Int
                      ( eval
                        ( U(Y),
                          Ans[j]
                        )(Y),
                        Y = 0 .. 1
                      )
                    );
      Q_b[j] := evalf
                ( Int
                  ( eval
                    ( U(Y),
                      Ans[j]
                    )(Y),
                    Y = 0 .. 1
                  )
                );
  end do;

  for j to numelems(MVals) do
      Ans1[j] := dsolve
                 ( eval([OdeSys1, Cond1], M = MVals[j]),
                   numeric,
                   output = listprocedure
                 );
      Theta_B[j] := evalf
                    ( Int
                      ( eval
                        ( V(Y)*Phi(Y),
                          Ans1[j]
                        )(Y),
                        Y = 0 .. 1
                      )
                      /
                      Int
                      ( eval
                        ( V(Y),
                          Ans1[j]
                        )(Y),
                        Y = 0 .. 1
                      )
                    );
      Q_B[j] := evalf
                ( Int
                  ( eval
                    ( V(Y),
                      Ans1[j]
                    )(Y),
                    Y = 0 .. 1
                  )
                );
  end do;