1) use alias

restart:
alias(pode = plots:-odeplot):
p:= dsolve({D(y)(x) = y(x), y(0)=1}, type=numeric, range=-5..2):
pode(p);

2) use macro

restart
macro(pode = plots:-odeplot):
p:= dsolve({D(y)(x) = y(x), y(0)=1}, type=numeric, range=-5..2):
pode(p);

Wait for more skilled people to explain you the pros and cons of alias versus macro

I don't understand why "normality" as something to do with that?
Let A the data (a list or a vector).

This gives the interval of variation and the range

(min..max)(A);
(max-min)(A);

Instead of the last command you can also use this one

Statistics:-Range(A)

If you want more specific things like inter quantile range or quantile ranges, look to these terms in the Statistics package.
InterQuantileRange, Quantile, Range are not dependent on the underlying statistical distribution, so normality doesn't matter.

I seriously doubt that your Analytical_sol is the solution of [ode, bc].

Here is the correct formal solution:
(LargeExpressions is used ot ease Maple to calculate the solution)

with(LargeExpressions):
isolate(ode, diff(theta(y), y$2)):
rhs(%):
compact := collect(%, y, Veil[C]);
CS := [ seq( C[i] = Unveil[C]( C[i] ), i=1..LastUsed[C] ) ]:

                12                   11   3779136       10
 68024448 C[1] y   + 136048896 C[2] y   - ------- C[3] y  
                                             5            

      1259712       9   157464       8   209952       7
    - ------- C[4] y  + ------ C[5] y  + ------ C[6] y 
         5                25               25          

      5832       6   5832       5    81        4   27         3
    - ---- C[7] y  - ---- C[8] y  + ---- C[9] y  + --- C[10] y 
      125            125            1250           125         

       9         2   3            1       
    + --- C[11] y  + -- C[12] y + -- C[13]
      100            20           32      


dsolve([diff(theta(y), y$2)=compact, bc]):
Compact_Exact_Solution := rhs(%): 

# The solution is a 14 th degree polynomial in y degree(Compact_Exact_Solution, y) 
# while your Analytical_sol is a 4th degree polynomial in y

degree(Compact_Exact_Solution, y);

                               14

I would like to point out that if you had read the answer I had given to your initial question (which you destroyed without informing me), you would have had the (correct) answer to this one too.

To_Read_before_throwing_away.mw

linestyle is not correctly defined, proceed as you do for color :
 

colors := [black, red, blue, green, orange, cyan, purple, yellow, navy]; 
styles := [solid, longdash, spacedash, dash, dashdot]; 

for j to numelems(k__vals) do 
sol := dsolve(
         eval(
           [ode, bc], 
           [k = k__vals[j], We = .1, Br = .3, x = 0, Q = Q__vals[j], lambda = lambda__vals[j]]
         ), 
         numeric
       ); 
plt[j] := plots:-odeplot(
            sol
            , [y, theta(y)]
​​​​​​​            , y = 0 .. 1 
​​​​​​​            , axes = boxed 
            , color = colors[j] 
          # , style = [line] 
            , linestyle = styles[j] 
          # , symbol = [solidcircle, asterisk] # useless when style=line
          ) 
end do; 
plots:-display([seq(plt[j], j = 1 .. 5)])

By the way, the ode can be integrated formally.
(with a little help to Maple's 2015 version)
Help_graph_style_mmcdara.mw

Customize theattached file as you want (font, line thickness, characteristic vertical values and so on)

 

Download pot_mmcdara.mw

From  https://mathworld.wolfram.com/SimpleGraph.html :
"A simple graph, also called a strict graph (Tutte 1998, p. 2), is an unweighted, undirected graph containing no graph loops or multiple edges (Gibbons 1985, p. 2; West 2000, p. 2; Bronshtein and Semendyayev 2004, p. 346). A simple graph may be either connected or disconnected.".

Refering to Tutte the notion of simple graph only applies to undirected graph.

In case you are not satisfied with RandomDigraph (like @dharr  I use Maple 2015 too) you can replace this command with these lines

restart:
with(LinearAlgebra):
with(GraphTheory):
A := RandomMatrix(10, 10, generator=0..1, density=0.2);
A := A - DiagonalMatrix(Diagonal(A)):  # set A[i, i]=0 to avoid loops

G := Graph(A):
DrawGraph(G);

In case you want to create an random undirected graph to match Tutte's definition,  add these lines:

G_undirected := Graph(Vertices(G), convert~(Edges(G), set)):
DrawGraph(G_undirected);

Itis very likely that the image you display is a 2D projection of a 3D surface.

Here is an example (I arbitrary choosed Q as the second independent variable, be free to select abnother one if you want to).

plot3d(
  subs(We = .1, n = .5, lambda = 0.95e-2, x = 0, k = .1, u) 
  , y = 0 .. 1 
  , Q = 0 .. .6 
​​​​​​​  , axes = boxed 
​​​​​​​  , labels = ["y", "Q", "u"]
​​​​​​​  , style = surfacewireframe
  , shading=zhue
​​​​​​​  , orientation = [0, 90, 90]
​​​​​​​)

Putting aside the issue which concerns the execution time...

Maple's DataFrame is not that far from R's dataframe. Maple's DataFrame has less capabilities than R's dataframe, but the "logics of manipulation" are quite similar. Of course syntaxes differ, but if you are familiar with R, using dataframes in Maple should be relatively simple.
Just keep in mind that not all R's capabilities are avaliable in Maple's DataFrame.

A main issue in using xls/xlsx files is that ExcelTools:-Import and ExcelTools:-Export are operations which are rather slow.
My advice is to use binary file as soon as the number of data is large.

For instance

series(F, t = infinity) assuming theta > 0, mu > 1

It's not impossible that series eturn a result with other assumptions about alpha and beta (?)

Let E some event whose the true but unknown probability of occurrence is noted p.
A classical problem is the assessment of p.

Frequentist way
One realizes N trials of the phenomenon which can make E to appear.
Let K the number of times E has been observed.
The (basic) frequentist estimation of p is f=K/N.

Bayesian way
One still realizes N trials and observe K occurrences of E.
Obviously the value of K is a realization of some Binomial random variable Q with parameters N and q (q being any value between 0 and 1). It's clear that some values of q are more likely than others.
The bayesian way relies upon the knowledge of the values Q might take, and represents this knowledge in termes of a specific statistic distribution named Prior.
For Objective Bayesians, the classical choice is Prior(Q) = Beta(a, b) where a and b have, at least in the non hierarhical bayesian approach, fixed values.

One can show that a Beta distribution is a conjugate prior of a Binomial likelihood, which means that the posterior distribution of Q belongs to the same family than the prior.
Thus the main result : Post(Q) = Beta(a+K, b+N-k)

See for instance https://en.wikipedia.org/wiki/Conjugate_prior

So there is no need to use a Metropolis Algorithm here.
With  a=21, b=1, N=4, K=2, Post(Q) = Beta(23, 3)

restart:
UseHardwareFloats := false
                             false
with(Statistics):
Mean('Beta'(23, 3));
evalf(%);
                               23
                               --
                               26
                          0.8846153846
Quantile('Beta'(23, 3), 0.025, numeric); 
                          0.7396941579
Quantile('Beta'(23, 3), 0.975, numeric);
                          0.9745346034

Remark: The convergence of your Markov Chain is slow because the your proposal (a Uniform(0, 1)) is highly inadequate.

I don't have time right now to go further but I bet that the posterior can be constructed in a thousand steps if the proposal is correctly chosen (as your prior concentrates in the range 0.8..1, take for instance a Uniform(-0.2/4, +0.2/4)). 

WATCHOUT : The prior being bounded, it is sure that some "prop" values you generate do not belong to [0, 1] (more certainly exceed 1). The expectation of the prior is 21/22 ~ 0.9545454545, meaning there are roughly 94% of chances that prop exceeds 1... look to this file
Download Beta+Uniform.mw

When you construct a Markov Chain an important quastion is "Did my chain behave correctly?"
A good looking chain must explore the space of possible values, must not be trapped on the same value for too many steps, its sample (after the burning phase) must posses a "correct" autocorrelation funcion and so on.
Just a question: have you observed the behaviour of your chain.

First: as you do not know a priori the number of iterations, it will be diffficult to predefine the array you talk about.
Secondly: the number of iterations depend on the method you use and on its parameters.
Thirdly: some methods have adaptive x steps and there is no simple relation between the number of the iteration and the physical time.

If yf1 and f2 do depend on the current iteration, why not write f1(x) and f2(x) instead ?

About events, look my answer to your previous question.

"Basically there is a very simply system of two differential equations and I would simply like to see an event triggered when the dependent variable is above 2"

In your worksheet you have written 

events = [[[0, 2-x < 0], [i(x) = i(x)+1]]]

but x is the independent variable.

Here are a few examples: the first and the second one correspond to what I understood (I arbitrarily used y(x) as the dependent variable used to monitor i(x)).

The third example is more complex: i(x) is augmented or reduce by one each time y(x) takes values +2 or -2. This example is purely illustrative.

The last example is a workaround to the condition "x > 2": it uses an auxuliary dependent variable w(x) defined by 

diff(w(x), x)=1, w(0)=0

Then, obviously, w(x)=x: "your" condition 

events = [[[0, 2-x < 0], [i(x) = i(x)+1]]]

then writes

events = [ [w(x)=2, i(x) = i(x)+1] ]

dsolve_events_mmcdara.mw

foo := proc(a, b) 
local x, y; 
x := a^2 - b; 
y := b^2 - a; 
printf("x=%a, y=%a, tests on return\n", x, y):
# or, if you prefer
# printf("%a, %a, tests on return\n", x, y):
return x, y
end proc:

A, B := foo(6, 3):
x=33, y=3, tests on return

A, B
                             33, 3

The simple way is 

ha := Histogram(
       samp
       , bincount = NumBins
       , size = [600, 400]
       , frequencyscale=absolute
       , gridlines=true
     ):
ha;

You can customise the plot this way

hb := plots:-display(
        ha
        , axis[2]=[tickmarks=[seq(i = nprintf("%2.1f%%", i), i = 0..8, 0.5)], color=red]
        , transparency=0.3
      ):
hb;

Histogram_percentage_mmcdara.m

Last improvement : in hb I set the max value of the y-axis to 8.
You can make this operation programatically:

pc_max := ceil(max(seq(plottools:-getdata(ha)[k][3][4, 2], k=1..NumBins))):

hb := plots:-display(
        ha
        , axis[2]=[tickmarks=[seq(i = nprintf("%2.0f%%", i), i = 0..pc_max, 1)], color=red]
        , transparency=0.3
      ):
hb;

... there is no way to do that with "standard" legend options.

Which doesn't mean it can't be done otherwise.
Here is a workaround: get inspired (if you want) and adapt it as you wish (if uou wtill want).

Remark: the values 0.95, 0.8, 1.2 and 0.03 have to be adjusted by hand

Download User_Legend.mw