Be carefull, this method doesn't handle correctly all the situations.

restart
alias (U=U(xi)):
Balance := proc(EQ)
  local eq   := eval(EQ, Diff=diff);
  local HODD := degree(eval(eval(eq, U=exp(__k*xi)), xi=0), __k);
  local HNLD := degree(eval(eq, U=__C), __C);
  M = solve(HNLD*M = HODD+M)
end proc:

eq := a*U+b*U^2+diff(U, xi$2)+diff(U, xi):
Balance(eq)
                             M = 2

eq := a*U+b*U^2+Diff(U, xi$2)+Diff(U, xi):
Balance(eq)

                             M = 2

eq := a*U+b*U^(-2)+diff(U, xi$5)+diff(U, xi$3)+c*U^3:
Balance(eq);
                                 5
                             M = -
                                 2

Let 

f := a -> x^3+a*x+1

a third degree polynom wrt x, parameterized by some real a.
In this simple case solve(f(a), x) doesn't return a list, nor a set, but a sequence:

f := a -> x^3+a*x+1solve(f(a), x):
whattype(%);
                            exprseq

If one defines a procedure sol this way

sol := a -> {solve(f(a), x)}:

one can easily verify that repeted executions of sol(a), where a is now assigned a numeric value, alwys returns theroots in the same order (defined by set):

1. If the three roots are real the order is defined by '<', that is the smallest root is the first oneand largest the last one.
2. If only one root is real and the two others are of the form p-q*I and  p+q*I (with q > 0), then the order is:
   1. the real root is the first,
   2. root  p-q*I is the second one,
   3. root  p+q*I is the last (third) one.

To avoid mixing the roots define the variable sola  this way:

sola := sol(a):

Then sola contains roots sorted according to the order given in points 1, 2, 3. This is because the formal result sol(a) contains a real root and two a priori conjugate complex roots (whoch might happento be both real for some values of a).

If you make a to vary in some real range, the first root in sola, which is assumed to be real, will take indeed real values, but also complex values, which means that its rank can take values 1, 2 or 3 depending on the value of a.
The plots in the attached file show how the rank of the three roots evolve as a ranges from -5 to +5.

To answer your question:

• The order of roots depends on their expression (real, p-q*I, p+q*I) and then a root who was ranked first for some value of a, can be ranked 2 or 3 for a close value of a, giving the false impression of a "random" order returned by solve.
• To understand (and control it?) this phenomenon you must "stuck" the roots, for example by doing sol := a -> {solve(f(a), x)}, and assign the roots a name, for instance

  R1, R2, R3 := sola[]:
  
  # check
  R1;
  R2;
  R3;

  Doing this will help you to track a given root as a evolves.

root_order.mw

Download Counter.mw

An alternative: search on this site for "progression bar" 

I suggest you to solve the pde without any boundary condition, which will givea solution that contains integration constants, and next to evaluate these constants given the boundary condition.

If I'm not mistaken thesolution should be this one: laplacian_mmcdara.mw

Here is an example of compartmental model (4 compartments) and a way to assess its transition coefficients given some observations.
As I do not have any observation for this hypothetical model I used the classical technique wich consists in creating "pseudo-observations" which correspond to the value of this model at some times.
More precisely for a N-compartments model M(x₁(t), ..., x_N(t), t ; K) with transition matrix K  (K numeric), the "pseudo-observations" are mape of P (N+1)-tuples { (t, x[1](t), ..., x[N](t)), t in {t₁, ..., t_P} }.
The goal is then to recover blindly estimations of K by exploting some prior domain Dom(K).
Without observation noise this estimation should be equal to K provided the model is identifiable (which depends on K and of the observations).

The attached file presents the case of noise free observations and a noisy case.

Hope this will give you some idea to go further.
An_example.mw
An_example_edited.mw

For a 9-compartment model you have potentially 90 coefficients to assess (9^2 compartment-to-compartment coefficients plus 9 compartment-to-exterior coefficients) which is almost impossible with usual methods unless the 9-by-9 matrix K is very parse (by "usual methods" I mean optimization methods by contrast to neural network methods for instance)
Even with specific estimation methods it's very unlikely that the solution you get will be "numerically unique": it often happens that the estimation which minimizes some discrepancy between model predictions and observations is not a point but a variety of higher dimension or, even if this is indeed a single point, that large deviations in some directions around it (think to this point as a point in dimension 90 for instance) do not change significantly the value if the discrepancy.

So be extremely carefull in applying a minimization algorithm.
A preliminary analysis may help. For instance, for large times, the behavior of some x_n(t) is generally governed by the largest coefficient k_{n, m}: this is a trick which can enable assessing the values of the leading coefficients of each equation?

Without knowing your model it's not possible to say more (see @Carl Love's reply).

I assume the histograms you refer two are those created by the package Statistics?
If they come from package ImageTools what I say doesn't apply.

Download HowToCatchTheBins.mw

In this kind of problem  the causality principle must hold: the Delta pulse is the cause of y(t), meaning that y(t) doesn't exist before Delta is applied.
So it doesn't matter to consider the limit of y(t) when t goes to 0 from the left.
The Laplace transform, for instance, is perfectly suited tosolve this (linear) problem:

restart
de := D(y)(t) + y(t)/tau = Dirac(t)/tau;
                             y(t)   Dirac(t)
                   D(y)(t) + ---- = --------
                             tau      tau   
with(inttrans):
lde := laplace(de, t, p)
                                   laplace(y(t), t, p)    1 
    p laplace(y(t), t, p) - y(0) + ------------------- = ---
                                           tau           tau
lde0 := eval(lde, y(0)=0):
ly := rhs(isolate(lde0, laplace(y(t), t, p))):
Y := unapply(invlaplace(ly, p, t), t)
        /   t \
     exp|- ---|
        \  tau/
t -> ----------
        tau    

Download laplace.mw

Two solutions:

restart:
with(LinearAlgebra):

sanity := proc()
  local A, __U, __S, __Vt;
  A :=RandomMatrix(3,10);
  __U, __S, __Vt := SingularValues(A, output=['U','S','Vt']);
end proc:

sanity():

or (IMO better for future use)

restart:
with(LinearAlgebra):

sanity := proc()
  local A;
  A :=RandomMatrix(3,10);
  SingularValues(A, output=['U','S','Vt']);
end proc:

U, VS, Vt := sanity():

The (A?) key to get a solution up to xi=0.25 is to use a continuation method (see here help(dsolve, numeric_bvp, advanced) ).
I wasn't capable to get a solution for higher values of xi.

Does this "critical" value of 0.25 means something to you?
If not you can try, as suggested by the error messages, to change the continuation strategy... but this likely require a lot of time that I have not to help you further.

secod_code_mmcdara.mw

in exponentials:

Download Num_PDE_mmcdara.mw

I agree with Rouben.
You can easily find (on the web or in your math books) the equation of the tangent plane to a sphere and try to code them in Maple.
Maybe this is what you have been asked for?

If you are not comfortable with maths and need to solve this problem quickly, you can use the geom3d package:
With_geom3d.mw
This will answer your question, but it won't help you understand how the equations were constructed.

Maybe reading carefully the three procedures TangentPlane, IsOnObject and onobjl could help you understand how to derive the equations.

here is an extremely simple way do draw random triangles
 

Download RandomTriangles.mw

Let me know if you have some constraints on the triangles.
For instance here is a simple way to generate triangles which contain the center of the circle:

 

Download RandomTriangle_2.mw

Note that the procedure in RandomTriangle_2.mw doesn't produce uniformly distributed triangles in the following sense: while p[2] is obviously uniformly distributed over [0, Pi], one can see that q (the third point p[3]) is not uniformly distributed over [Pi, 2*Pi]

# Do this to convince you
K := 10000:
P := Matrix(K, 2):
for k from 1 to K do
  p := rand(0. .. 1.*Pi)():
  q := rand(1.*Pi.. p+1.*Pi)():
  P[k, 1] := p:
  P[k, 2] := q:
end do:
Histogram(P[..,1]);
Histogram(P[..,2]);

Here is the exact expression of the PDF of q:

p2 := RandomVariable(Uniform(0, Pi));
q  := Pi + RandomVariable(Uniform(0, 1)) * p2;

# and thus 
combine(PDF(Y, y), ln);
             /                                 /  -y + Pi\  
             |                               ln|- -------|  
             |                                 \    Pi   /  
    piecewise|y - Pi < 0, 0, y - Pi <= Pi, - -------------, 
             \                                    Pi        

                    \
                    |
                    |
      Pi < y - Pi, 0|
                    /

(it's likely that using the remember option would be more astute)

I used a predator velocity equal to the prey velocivy divided by the eccentricity of the ellipses (see me reply under the login sand15)

Download Capture.mw

to get he outputs you want.
Among them:

Download Riemann_sum.mw

A help page that might interest you

help(Student:-Calculus1:-RiemannSum);

Riemann_sum_2.mw

Finally, better presentations (IMO) can be obtained, for instance with DocumentTools:-Layout.

Drawing a parabolic cylinder with generatrix paralel to axis y (for instance) is easy:

restart
with(plots):

implicitplot3d(
  y-x^2
  , x=-1..1, y=0..1, z=0..1, grid=[40$3]
  , style=surface, color=cyan, transparency=0.3
  , axes=normal
  , seq(axis[k]=[tickmarks=0], k=1..3)
);

Drawing the magenta curve requires you define what this curve is.
Once done: unfold the parabolic cylinder (this gives the plane y=0) and draw the curve C on this plane; next "fold" this plane to retrieve the original parabolic cylinder.
The curve C' you get is the image of C, drawn on the parabolic cylinder.
The parametric equation of curve C' is

[x, f(x), g(x)] for any x in some drawing interval

Where:

• f := x -> a*x^2 + b*x + c represents the equation of any cross section oc the paramolic cylinder;
• z = g(x) is the explicit equation of curve C in plane y=0.

Here is an example when the curve C is given by g : x -> g(x)=x
pc0.mw

The file above contains another example

pc.mw