First solution : click on the hand right to the exclamation point (it should be red)

Sometimes this doesn't work if the build in function (here fsolve) has run some child process.
In this case the only thing to do is wait, hopping the command will stop correctly ... or to kill the process itself.

Second solution: use kernelopts(cputime=.....).
For instance kernelopts(cputime=15) will stop your active worksheet as soon as the cpu time reaches 15s.
Watch out: kernelopts(cputime=15) is global to the worksheet: if 10 s has already been used for some computations, the remaining ones have to be executed in less than 5 s  to complete entirely.
When cpulimit is reached (I cannot guarantee it'd the same thing on all platforms or systems), you may receive a "kernel connection lost" message. If it so, you will not be able to do any other thing within your worksheet (other worksheets are still operational) and the only thing tou will be able to do is to save the active worksheet, to close it and to open it again.

cpulimit.mw

I don't know if Maple offers a softer mechanism to kill a command programatically. 
 

The first part of the attached file is essentially a corrected version of yours
In the second part you will find some improvements that could help you to do the same things in a simpler way

Could_help.mw

Download YetAnotherAlternative.mw

If I concentrate on your sentence "I am trying to generate a binomial sample", and if you want to use the student package, then the solution is obvious:


with(Student[Statistics]):
N := 50:
P := 0.1:
X := BinomialRandomVariable(N,P):

# S is a sample of size M 
M := 10:
S := Sample(X, M);

But maybe, as Carl Love sugegsted, you want to do something different ?


PS : for my own I prefer to use the Statistics package.
In this case :

restart:
with(Statistics):
N := 50:
P := 0.1:
X := RandomVariable(Binomial(N, P)):
M := 10:
S := Sample(X, M);

WATCH OUT: here S is a vector of floats (not integers as with Student[Statistics]): be aware of this when you use S! 
 

restart:
M := 5;
A := Matrix(M+1, M+1, (r,s) -> `if`(type(r+s, odd), 2^(k+1)*sqrt((2*r-1)*(2*s-1)), 0))

restart:
R := ((n1-n2)/(n1+n2))^2:
simplify(subs(n1=xi*n2, R));
plot(%, xi=-1..2, labels=[n1/n2, 'R']);

If you want a symbolic solution just use
dsolve({ode, bcs}, f(x))

But your second bc is not correctly stated : you should write f((4*10^(-8)-rho)^.5)= "some value"

Unfortunately, dsolve will return no result because the ode is not tractable symbolically.

Then you can try dsolve({ode, bcs}, numeric).
In which case you need to fix initial conditions, not "noundary" ones.
More of this, you need to fix the value of rho (unless you use paramerized dsolve : dsolve({ode, bcs}, numeric, parameters=[rho]) )

the correct syntax is a:=a+1  not (a=a+1) (idem for b, c, d)

answer.mw

First remark :
I loaded your file  Advheat-53.mw and got this error
Error, recursive assignment

It comes from the line "eqs := eqs, ....  " and tells you that "eqs" has not been already defined.


Second remark
In Advheat-54.mw : solve operates on a single equation or a set of equation (at least until MAPLE 2016 ... I do not have more recent versions).
The syntax of your last line is not correct and should be 
Array(1 .. 9, rhs~(convert( solve(eqs, {seq(T[k], k = 1 .. 9)}) , list)))

In effect, eqs has a lot of indeterminates and you must tell MAPLE what are the ones it has to solve for, this is the meaning of {seq(T[k], k = 1 .. 9)}.

Now a computational issue: even if you execute this command you will be very disappointed by the result you'll get. I'm not even sur that you will obtain a solution (I can't).
Do not forget you have several fourth order equations to solve: the result will be amazingly huge if you look for a formal one.
To convince you, replace all the "4" in the fourth power by "m" in your eqyuations and execute 

Array(1 .. 9, rhs~(convert( solve( subs(m=1, eqs), {seq(T[k], k = 1 .. 9)}) , list)))

(ok, this is not the correct exponent, but just look to the result)

Now do the same thing for m=2 and observe what happens.


Some hints

1. If a numerical solution satisfies you, I suggest to fix the values of all the quantities, but the T[1..9], in eqs and use fsolve instead of solve.
   Syntax    Array(1 .. 9, rhs~(convert( solve(eqs) , list)))
    
2. If you need want a formal solution of this advection problem, MAPLE proposes you the procedure "pdsolve". Please have a look to it and let me know if it could suit you.
    
3. It's likely that someone here will have other(and better) ideas

But for the moment I'm affraid that yout problem, in the way you stated it, is untractable (?)

And a third one :

1 +~ Statistics:-CumulativeSum(X)

You see here that MAPLE has many overlapping feature: look to Carl Love's answer.

BTW : I didn't know the procedure PartialSums of the ListTools package

@Rouben Rostamian  @acer  @Mariusz Iwaniuk 

Some of the solutions are based on some arbitrary cutoff and I do not consider they are correct.

Rouben's cutoff-free solution is far better ... if we accept a zero gravity of course !!!

You will find in the attached file another solution free of any trick  (excepted the instruction in pink used here to prevent null time ranges).
The solution is piecewise constructed (a strategy I use to use for more complex situations than the bouncing ball test case).

Make a zoom around the 16th rebound (t=4.2425543703...) for a detailed comparison between Rouben's solution and mine's


bouncing_ball.mw

 

 1/ There is no deed to write f(x_i,y_i)=0.01 explicitely: just define a new function g(x_i, y_i) = f(x_i, y_i)-0.01
      The "constraint" the becomes g(x_i, y_i) = 0. This is not simpler: we just discard a useless numerical value.

2/ You say that f depends on some parameter, let's say p, you "calculated" by some optization process. 
     When does this process occur? Befor generatind the points, or after having generating them ?
     
Without more explanation I'm not sure anyone can give you the answer you are looking for?


Example   :  g(x, y, p) = x-y-p
                    Your "constraint" is g(x, y, p)=0.
                    A trivial solution is x_i = y_i for all i and p = 0
                    Another one is x_i = y_1 + 1 and p = -1

Please have a look to the attached file.
I do not find your second-order approximation but some mistake of my side is not exluded (it would have been safer for me to provide Maple code instead of Latex code).

qwerty.mw

First point : Error, (in int) integration range or variable must be provided

pn:=proc(i,n,t)
if n=1 then return int(hd[i],t) fi:
return int(pn(i,n-1,t)):
end proc:

When n = 2 (your instruction  tmp := alpha1(t)*pn(i,1,t)+alpha2(t)*pn(i,2,t): ) the error message means there is no integration variable.
Changing int(pn(i,n-1,t)) by int(pn(i,n-1,t), t) eludes this message.


Second point: the execution of the line  A := Transpose(LinearSolve(Transpose(H+TMP), R))  takes a very huge amount of time (Maple 2017.2).
So I propose to rewrite the system in a more symbolic way (see the attached file   # let's try to solve the system in a more symbolic way)

Third point   Error, missing operator or `;`
Change "sum" into "add" after line #Now compute the approximate solution

example_1.mw
 

@raghav6594 

Without presuming about what you really want to do ...

Your function is zero almost everywhere, which will be very difficult, if not impossible, to use it in an ODE/PDE problem.
Maybe a softened version of it could be sufficient ?
Without any pretention the attached file presents a simple idea to implement and use such a softened function.

sha.mw