Either of these work in Maple 2015

  restart;
  WM1 :=proc( P, n)
             plots:-matrixplot(LinearAlgebra:-MatrixPower(P, n), heights=histogram):
        end proc:
  U := Matrix(2$2, [0.8, 0.2, 0.4, 0.6]):
  plots:-display( seq(WM1(U,n), n=1..10), insequence=true);

or

 restart;
  WM2:=proc( P, n)
              LinearAlgebra:-MatrixPower(P, n):
       end proc:
  U := Matrix(2$2, [0.8, 0.2, 0.4, 0.6]):
  WM2(U, 1);
  plots:-animate(plots:-matrixplot, [WM2(U,n), heights=histogram], n=[$1..10]);

NB the use of the LinearAlgebra:-MatrixPower(P, n) is not required but (for matrices) is generally a lot more efficient than using P^n

well the attached (expPrd.mw) works in Maple 18 ( ie released in 2014) and allsubsequent Maple releases. I'm afraid I can't verify in any older Maple releases: I only keep about the last eight (or so!)  'major' Maple releases running at any one time, and Maple 17 has "dropped of the list"

Implemented code is

  restart;
  kernelopts(version);
  a := {x+1, x+2, x^2+1, x^2+x+2}:
  seq([expand(`*`(j[]))], j in  combinat:-choose(a, 3));

You have a few problems

1. As far as I can tell, only numerical solutiions are possible so numeric values for 'k' and 'infinity'. In the attached ( quadBVP.mw ), I have used k=1 and infinity=10 (but also see comments in (4)  below)
2. The highest derivative in the second ode of your system is quadratic. This mean that the system cannot be unambiguously converted to a first-order system (becuase the quadratic has two roots). One can compute the two possible roots analytically, substitute the selected root into the ode system, then use dsolve(..numeric)
3. For the parameter values selected in (1) above, one of the roots obtained in (2) above leads to real answers, and the other fails because complex numbers are produced. (It may be possible to split everythig into real and imaginary parts then solve this extended system, but I have made no attempt to do this in the attached.)
4. Out of idle curiosity, if I use the parameter values k=1, infinity=100, then both roots obtained at stage (2) above lead to complex solutions at step (3) above. On the other hand, if I use the parameter values k=1, infinity=1, both both roots obtained at stage (2) above lead to real (and different) solutions in step (3) above

The code shown

  restart;
#
# Solution will only be obtained using numeric methods.
# Need to decide numeric values for 'k' and "infinity"
# The following choices are arbitrary - be careful other
# choices *may* lead cuase the following solution process
# to fail
#
  k:=1:
  Inf:=10:
#
# Define the odes and the boundary conditions
#
  odes0 := [   3*diff(f(x), x, x, x)
             + 18*k*diff(f(x), x, x, x)*diff(f(x), x, x)^2
             + 2*diff(f(x), x, x)*f(x)
             - diff(f(x), x)^2
             = 0,
               3*diff(g(x), x, x)
             + 3*diff(g(x), x, x)^2
             + 2*diff(g(x), x)*f(x)
             - 2*diff(f(x),x)*g(x)
             = 0
          ]:
  bcs:= [ f(0) = 0,
          D(f)(0) = 1,
          D(f)(Inf)=0,
          g(0) = 1,
          g(Inf)=0
        ]:

#
# First (brain-dead) attempt at a solution. NB this will
# result in an error
#
  dsolve( [odes0, bcs[]], numeric);
#
# The reason for the above error is that the second ODE
# in the system is quadratic in the highest derivative
# diff(g(x),x,x). As a general rule quadratics have two
# solutions - which one to pick???
#
# Treating the second ode above as a simple quadratic in,
# diff(g(x),x,x), get an analytic solution for this
# quantity solve. Obviously two roots are returned 
#
  rts:=solve( odes0[2],
              diff(g(x), x, x)
            );
#
# Quite arbitrarily, select the first root obtained above
# and substitute this for diff(g(x), x, x) in the ODE
# sysstem
  odes1:=eval( odes0,
               diff(g(x),x,x)^2=rts[1]^2
             ):
#
# Now solve the system
#
  sol1:=dsolve( [ odes1[],
                  bcs[]
                ],
                numeric
              ):
#
# Plot the solution
#
  plots:-odeplot( sol1,
                  [ [x, f(x)],
                    [x, g(x)]
                  ],
                  x=0..Inf,
                  color=[red, blue],
                  title=typeset( "%1 and %2 vs  %3 for k= %4 and infinity= %5",
                                 f(x), g(x), x, k, Inf
                               ),
                  titlefont=[times, 16],
                  legend=[ typeset(f(x)),
                           typeset(g(x))
                         ]
                );
#
# An obvious question - what happens if one uses the
# other root of diff(g(x),x,x) obtained above??
#
# This produces an error because complex numbers are
# returned - so maybe the is the answer the OP doesn't
# want/need!
#
  odes2:=eval( odes0,
               diff(g(x),x,x)^2=rts[2]^2
             );
  sol2:=dsolve( [ odes2[],
                  bcs[]
                ],
                numeric
              );
  plots:-odeplot( sol2,
                  [ [x, f(x)],
                    [x, g(x)]
                  ],
                  x=0..Inf,
                  color=[red, blue],
                  title=typeset( "%1 and %2 vs  %3 for k= %4 and infinity= %5",
                                 f(x), g(x), x, k, Inf
                               ),
                  titlefont=[times, 16],
                  legend=[ typeset(f(x)),
                           typeset(g(x))
                         ]
                );

leads to the figure

many variations  on Black-Scholes option pricing in Maple's finance package, so for example the code

with(Finance):
plot( [ BlackScholesPrice(30, 30, t, 0.2, 0.1, 0, 'call'),
        BlackScholesPrice(30, 30, t, 0.2, 0.1, 0, 'put')
      ], 
      t=1..10,
      color=[red, blue]
    );

will produce the plot

and quite a nasty one, because I think we can all agree that

restart;
type(x^a*y,  `&*`(identical(x^a), name));
type(x^2*y,  `&*`(identical(x^2), name));
                              true

                             false

should return true, true. I will submit an SCR

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

I *think* (but obviously cannot prove)  that the origin of the problem lies with the Maple internal data structures PROD and POWER, and the "translation" between these and "user" representations.  (Please don't confuse the word ;structure' in the phrase 'internal data structure' with the word 'structure' associated with a Maple type - same word but different meanings in these two contexts!)

I don't believe that any Maple user should have to know anything about Maple's internal data representations in order to use the programme. In particular any "translation" between Maple's internal data representation and "user representation" ought to be bulletproof,  so the following is 'just for interest'

If one examines the description  of Maple's internal  data structures in Appendix 1 of the the Maple Programming Guide.

From the description of the PROD internal data structure, (emphasis added)

This structure is interpreted as pairs of factors and their numeric exponents. Rational or integer expressions to an integer power are expanded. If a rational constant is in the product, this constant is moved to the first entry by the simplifier. A simple power, such as a², is represented as a PROD structure. More complex powers involving non-numeric exponents are represented as POWER structures.

Similarly the description of the POWER internal data structure, (emphasis added)

This structure is used to represent a power when the exponent is not an integer, rational, or floating-point value. When the exponent is numeric, the POWER structure is converted to a length 3 PROD structure.

These "internal" data structures can be confirmed using the dismantle() command - as in the following. Note that simple powers such as x² are indeed converted to a PROD structure - the second example below

dismantle(x*y);

PROD(5)
   NAME(4): x
   INTPOS(2): 1
   NAME(4): y
   INTPOS(2): 1


dismantle(x^2);

PROD(3)
   NAME(4): x
   INTPOS(2): 2


dismantle( x^2*y );

PROD(5)
   NAME(4): x
   INTPOS(2): 2
   NAME(4): y
   INTPOS(2): 1


dismantle( x^a*y);

PROD(5)
   POWER(3)
      NAME(4): x
      NAME(4): a
   INTPOS(2): 1
   NAME(4): y
   INTPOS(2): 1

is to use eval(variableName, answerSetOrList) - as in the attached

restart:
 CD1 := -.5: CD2 := 2: CD3 := 1.: 
 g1 := 5.: g2 := 3.: g3 := 2.: 
 cg1 := 0.: cg2 := 3.7: cg3 := 1:
 cd1 := 7.5: cd2 := 0: cd3 := 0: 
 L1 := .72: L2 := 8.6: L3 := 5.5: L4 := 1.25: L5 := 3.102: L6 := 5.1: L7 := 7.: 
f1 := (CD1-x4)^2+(CD2-x5)^2+(CD3-x6)^2-L1^2: 
f2 := x1-5: 
f3 := x5-1.45: 
f4 := (g1-x1)^2+(g2-x2)^2+(g3-x3)^2-L3^2: 
f5 := (x7-x1)^2+(x8-x2)^2+(x9-x3)^2-L2^2: 
f6 := (cg1-x7)^2+(cg2-x8)^2-(cg3-x9)^2-L4^2: 
f7 := x7+x8-1.2*cg2: 
f8 := (x4-x7)^2+(x5-x8)^2+(x6-x9)^2-L5^2: 
f9 := (cd1-x10)^2+(cd2-x11)^2+(cd3-x12)^2-L6^2: 
f10 := (x1-x10)^2+(x2-x11)^2+(x3-x12)^2-L7^2:
f11 := x11-.1*x12: 
f12 := x1-x2-x3-x4+x5+x6+x7+x8+x9+x10+x11-x12: 
ans1:=fsolve({f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12});
ans2:=RootFinding[Isolate]([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12], [x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12])[4];
eval(x1, ans1);
eval(x1, ans2);

since both of these procedures contain syntax errors - then they are not correct, by any definition

The problem is associated with the use of punctuation characters (`) in print statements

The attached code ( badprocs.mw ) at least executes - I just commented out all the print sttements

Whether or not it produces the answers you desire, I have no idea, since I have no knowledge of what calculation ou are trying to perform. See the comments by mmcdara

If you change Digits within a procedure, the new setting will be applied when the procedure is executing, but this change is undone when the procedure exits

so the code

 Digits:=2;
  f := proc()
       Digits:=Digits + 1; 
       1.0/3;
       end proc:
  f();
  Digits;
  1.0/3;

sets Digits equal to 2 at the top-level. The command f() executes the procedure: within this execution Digits will be set to 3, 1.0/3 will be evaluated to 3 digits, but this change to the setting of Digits will be undone when the procedure has finished executing.

Hence the aove code will return

                          Digits := 2

                             0.333

                               2

                              0.33

A small warning. The comments made by vv are correct in that in the command

Digits:=3;
  evalf[8](1.0/3);

Maple's automatic simplification will 'simplify' 1.0/3 using the current setting of Digits (ie 3) before executing evalf[8]() - so you are only ever going to get 3 digits.

On the other hand. if you were to write

Digits:=8;
  evalf[3](1.0/3);

the automatic simplification will use 8 digits, giving 0.33333333, and eval[3]() will evaluate this to 3 digits - so you will get 0.333

I transferred my answer and then deleted

232573-Error-In-Parametric-Plot

one trivial to fix, and one difficult

The trivial problem is that you assign your ode system to the name 'ODESs', but when you refere to the ode system within the procedure getPlot(), yuo use the name 'ODEs' - which (of course) is a completely different (undefined) quantity.

If I fix this, you code runs, but the odeplot() command fails with the error message

Warning, cannot evaluate the solution further right of .63668160e-2, probably a singularity

In the attached (odeProb.mw) I have restricted the plotting range to .6e-02, just to examine what is happening to the solution near the point where Maple is complaining about a singularity. This produces the plot shown below

This is not conclusive, but does suggest that B(T) is heading for infinity as T gets a little more than 0.006. Discovering precisely why ypur ODE system returns what appear to be "infinite" answers for quite small values of T is more difficult. A few possibilities

1. Your ODE system is stable, but there is a typo somewhere in your equations
2. Your ODE system may be stable for certain values of the parameter list 'params' and unstable for other values. If you believe this to be the case then it may be worth trying a few other sets of values in the list 'params'
3. Your ODE system is unstable. If the ODE system is modelling some "real-world" physical phenomenon, then ths is possible but unlikely. Very few "real-world" problems lead to "infinite" answers, unless the 'model' is inadequate in some way

one trivial to fix, and one difficult

The trivial problem is that you assign your ode system to the name 'ODESs', but when you refere to the ode system within the procedure getPlot(), yuo use the name 'ODEs' - which (of course) is a completely different (undefined) quantity.

If I fix this, you code runs, but the odeplot() command fails with the error message

Warning, cannot evaluate the solution further right of .63668160e-2, probably a singularity

In the attached (odeProb.mw) I have restricted the plotting range to .6e-02, just to examine what is happening to the solution near the point where Maple is complaining about a singularity. This produces the plot shown below

This is not conclusive, but does suggest that B(T) is heading for infinity as T gets a little more than 0.006. Discovering precisely why ypur ODE system returns what appear to be "infinite" answers for quite small values of T is more difficult. A few possibilities

1. Your ODE system is stable, but there is a typo somewhere in your equations
2. Your ODE system may be stable for certain values of the parameter list 'params' and unstable for other values. If you believe this to be the case then it may be worth trying a few other sets of values in the list 'params'
3. Your ODE system is unstable. If the ODE system is modelling some "real-world" physical phenomenon, then ths is possible but unlikely. Very few "real-world" problems lead to "infinite" answers, unless the 'model' is inadequate in some way

1. You can convert the expression displayed in your original post to a polynomial of degree 18
2. In general no analytic (ie "exact") solution exists for any polynomial with degree>4
3. However fsolve() with appropriate options will return all 18 "roots" - some of these are complex, which apparently you want to discard.
4. Of the remaining real roots some have multiplicity>1
5. If the real roots are restircted to those which are >2 (not 100% sure if this is a requirement??), then there appear ot be two roots, 2.1213203436 with a multiplicity of 3 and 2.5579901741 with multiplicity 1
6. There is a second problem in the worksheet you supply - a fifth order polynomial, so again no analytic solution is possible, although fsolve() with appropriate options will return five roots of which only one is real ,6.5433137308

See the attached roots.mw whihc implements the code

  restart;
  Digits:=20:
  interface(displayprecision=10):
  c := -(2*t)/3:
  s := x -> -x^3 - t*x^2:
  s(c):
  s(s(c));
  expr:=s(s(c))=(-t+sqrt(t^2-4))/2;
  expr2:= expr-op([2,1],expr);
  expr3:=expand(expr2^2);
  expr4:=lhs(expr3)-rhs(expr3);
  degree(expr4);
#
# Get all the solutions including the complex ones
# Since the polynomial is of degree 18 - there ought
# to be 18 of these (although some may have
# multiplicity >1)
#
  sols:=[fsolve(expr4, complex)];
#
# Check the number of solutions - yes there are 18
#
  numelems(sols);
#
# Select only the real solutions
#
  rsols1:=select( type, sols, realcons);
#
# From the above, select only thos solutions which are >2
# (not sure if/why this criterion exists!)
#
  rsols2:=select(j->j>2, rsols1);
#
#
# As a double check evaluate the expression at
# the real solutions which are >2. All answers should
# be zero, within floating point accuracy limits
#
  seq( evalf(eval( expr4, t=j)), j in rsols2);
  restart;
  Digits:=20;
  interface(displayprecision=10):
  expr:=128*x^5 - 288*x^4 - 1296*x^3 - 11664*x^2 - 13122*x - 59049 = 0;
#
# Get all the solutions - there ought to be five
#
  sols1:=[fsolve(expr, complex)];
  numelems(sols1);
#
# Select only the real solutions
#
  rsols1:=select( type, sols1, realcons);

The code

restart
f:=x->2*2^x-2;
g:=x->-1/2*x^2+3/2+5;
sol1:=fsolve( 2*2^x-2 = -1/2*x^2+3/2+5, x);
sol2:=fsolve( 2*2^x-2 = -1/2*x^2+3/2+5, x, avoid={x=sol1});

will return

                      sol1 := 1.787210719

                      sol2 := -4.094614685

I'm not sure I'd call that a "piecewise" boundary condition, but one way to produce the solution 3.1 in your reference is shown in the attached diffusion.mw which uses the code

restart
PDE:= diff(f(x,t),t)=Db*diff(f(x,t),x$2);
IC:= f(x,0)=Q*Dirac(x);
sol:=pdsolve([PDE, IC]) assuming t>0, Db>0;
#
# Producing the *exact* form required needs some
# "simplification" commands - and there may be
# more elegant ways than the following
#
  solSimp:=lhs(sol)=numer(rhs(sol))/combine(denom(rhs(sol)), symbolic);

entries(x, pairs)

doesn't supply what you want, then yu are going to have to be a bit more specific about the form of output you require - it will (probably) be provided using appropriate options to either/both of the commands entries() indices()