It works well without any corrections with Maple 2015.2.

A tip: try to use the package "Statistics" instead of "Student:-Statistics".
Then change 

B := BetaRandomVariable(61, 2)

by

B := RandomVariable(BetaDistribution(61, 2))

@9009134 thanks

Tracking the site you mention (Andrews' paper) I found this one where Mathematica has been used to plot the streamline  potential-flow-around-a-torus-laplace-equation-in-toroidal-coordinates.
Are they the kind of plots you woulf like to obtain with Maple?

 

@nm 

Much simpler than the solution I've obtained, great.

I agree with @rlopez , the figure you present has nothing to do with the laminar flow of a viscous fluid around a cylinder.
See here for instance petitRe.pdf

Where do your formulas come from?
 

@Mac Dude 

I agree, I guess "0.78 kg CO2/ kg" must be interpreted as "0.78 Kg of CO2 per kilogram of <something>".
So its seems to be more a question of text formatting than a question of units.
Maybe something like this (with Maple 2015 ; could be inproved in latter versions by using neutral operators %* and %/ instead of `per kg`)

0.78*Unit('kg')*`CO__2`*`per kg`

@Carl Love @acer

Your two proposals are not strictly equivalentfor the depend on the structure of the data lists.
Acer's seems to be a little bit mode generic.

Here LL represents carl's list of the lists of "y" values and LLL the list of the lists of "(x,y)" couples.
 

Download DataList.mw

Your model is a  linear model under a proper reparameterization.

Note that, even if it  is a quite common way to linearize a model, the method can lead to numerical instabilities if the auxiliary variables of the reparameterization are of very different magnitudes.

More of this, even in its original rational form the model depends only on 4 parameters (you can divide both its numerator and its denominator without changing the results).
 

restart;
with(Statistics):
X:=Vector([0,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.012]):
Y:=Vector([1.103,1.057,1.016,0.978,0.94,0.91,0.88,0.85,0.826,0.8,0.778,0.735]):
# Transform the original Model into an Alternate one named AltModel

Model := y-(-a*t+b)/(-d+e*t+f*t^3)=0;
normal(Model);
AltModel := numer(lhs(Model)); # everywhere denom(lhs(Model)) <> 0
                 
# Rename z__p each product (t^p)*y
# Then the model is a linear model of the form y=A + B*t + C*z__1 + D*z__3
# It can be found using LinearFit once the vectors corresponding to z__1 and z__3 have been constructed

ZModel := [1, t, z__1, z__3]:
res := Statistics:-LinearFit([1, t, z__1, z__3], <X | X*~Y | (X^~3)*~Y | Y >, [t, z__1, z__3])

# Express this resultl in native variables t and y

RES := y = eval(res, [z__1=y*t, z__3=y*t^3])

# Rebuild the model under rational form

isolate(RES, y);
Fn := rhs(%):

plots:-display(plot(X,Y,style=point,color=red,symbolsize=12), plot(Fn,t=X[1]..X[-1],color=blue), labels=["X","Y"]);

Download LinearFit.mw

@acer 

Excellent!

Just a question, why does 
solve(convert~(Or(x::A,x::C, x::B),relation)); 
return x while
solve(x<+infinity)
returns 
RealRange(-infinity, Open(infinity)) ?

@Ivi 

The ides is the following :

• Set u(t)=0.007-f(t)
   
• Then the equation writes diff(u(t),t) = beta*C1*u(t)^2 + C2*u(t)^3 where C1 and C2 are constants and beta is a short for A(t) = (1.75*exp(0.06*t)-0.75*exp(-200*t)) and, at this point, beta is considered as a constant.
   
• Set v(t)=1/u(t)
   
• Derive the EDO in v(t)
   
• Set the initial condition under the form v(T)=K, where T >=0 and K an appropriate constant.
  For instance, if T=0, then K=1 / 0.007 given the two transformations f(t) --> u(t) --> v(t).
   
• Solve it exactly (not numerically) : this is possible because (with beta constant) the ode is of a classical form.
   
• Let "sol" this exact solution (it depends on t, beta, T, K)
   
• Now proceed this way to construct the global solution over [0, tau] where tau > 0:
  • vsol := [ [0, 1/0.007]]
  • for theta from 0 to tau by delta_theta do
        beta := A(theta);
        T      := theta;
        K      := vsol[-1][2]; # value of v at the last generated point
        compute sol(beta, T, K, t=theta+step)  ; # means "advance" the solution from theta to theta+step
       add the point [theta, sol(beta, T, K, t=theta+step)] to vsol
    end do
     
• look vsol : the values of v are complex but the first ones have a very small imaginary part one can think is due to numerical precision.
  Keep only the points [t, v(t)] for wich Im(v(t)) is small.
   
• Plot the result. It appears to be very close to the result obtained with rkf45 for instance.
  This prooves the singularity diagnosed by dsolve[numeric] is not some numeric artefact, but a consequence of the LambertW function that the exact "local" (meaning "with constant beta") solution contains.

Download modest_contribution_2.mw

@Ivi 
I agree with Preben, 
This is what I used to do when Preben was working on his side. I had "identified" the coefficient befor diff(f(t),t) as the "critical" one.

 

Download modest_contribution.mw

@dharr 

I guess you refer to the sentence "I have reason to expect that this pattern will repeat for successive points (i.e. points 13, 14, ..., 30)"?
I wasn't sure to understand it correctly, in particular that [13, y[13]] was equal to [1, y[1]].

This is why I said that the OP would have specified either   [ [0, p[-1][2]], p[ ] ]   or  [ p[ ], [13, p[1][2]] ]
 

@erik10

Time and memory comparisons
 

restart
f1 := proc(Bins, Probs, N)
  local X, S;
  uses Statistics:
  UseHardwareFloats := false:
  X := RandomVariable(ProbabilityTable(Probs));
  S := round~(Sample(X, Vector(N, datatype=float))):
  Bins[convert(S, list)]
end proc:

f2 := proc(Bins, Probs, N)
  local X, S;
  uses Statistics:
  UseHardwareFloats := false:
  X := RandomVariable(EmpiricalDistribution(Bins, 'probabilities'=Probs)):
  S := Sample(X, N);
  round~(S)
end proc:

Probs := [0.10, 0.30, 0.20, 0.25, 0.15]:
Bins  := <50, 60, 70, 80, 90>:

S1 := CodeTools:-Usage(f1(Bins, Probs, 10^4)):
S2 := CodeTools:-Usage(f2(Bins, Probs, 10^4)):
memory used=13.89MiB, alloc change=32.00MiB, cpu time=120.00ms, real time=120.00ms, gc time=13.84ms
memory used=18.76MiB, alloc change=0 bytes, cpu time=96.00ms, real time=96.00ms, gc time=0ns
Statistics:-Tally(S1);
Statistics:-Tally(S2)
    [50 = 1005, 60 = 2989, 70 = 2035, 90 = 1428, 80 = 2543]
    [50 = 1012, 60 = 3033, 70 = 2076, 90 = 1442, 80 = 2437]

Download Comparison.mw

@Carl Love @vv  @Axel Vogt

Maybe not as simple as that for the original paper seems to contain a lot of typos (a lot of them come probably from digitalization of the original paper).
For instance, in procedure "av" the line 
for i: = i— 1while co x eps2 < p A i > 0 do
co x eps must probably be intended as co * eps  and  A  as /\ (for "and").
But what about for i: = i— 1whose syntax is not correct (even if "—" one must likely read as "-")

Given the problem (multimodal optimization) I understand "search of the global minimum (minima?) of a function with several local minima". For Axel: I'm not familiar with the "DirectSearch" package, does-it handles this kind of problem?
Ideally, but this is non free toolbox, GlobalOptimization seems more suited.

The fact the function is noisy was an anusual point by the time the original article has been published (1980). 
But it is today a subject widely investigated, and this  for at least 20 years, and a lot of efficient methods do exist to address this problem (see for instance scholar by the same author antanas žilinskas).
So I think it's probably better and simpler to implement in Maple a "modern" algorithm (or to use and adapt something already existing in Maple) instead of this old algorithm that we are not even sure it doesn't contain more ambiguous typos.

@dohashi  @Carl Love

Parallel programming being recurrently mentioned here, don't you think that all the courses, tips, advices, answers that both of you, and a few others, have done on the subject of parallel programming would be worth compiling?
Mapleprime's advanced search can help find useful information, but, with the exception of Dohashi's courses, it remains (IMO) too scattered to be really effective.

Pretty work!

When I was studying mechanics, years ago, I discovered Timoshenko's books and I still have the one I bought for a research project.
It was entitled "Théorie des vibrations" (in french, translates as "Vibration Theory"), and your work prompted me to reopen it. There's a whole chapter on "Lateral Vibration of Prismatic Beams". At that time (1954), the term "Krylov-Duncan" did not seem to be used (at least Timoshenko did not use it in this book).
Your work takes me back years and I read it with great pleasure.