You define f[i] as a function of two arguments p2 and pr but you use only one in the following command:

Cs := {seq(beta[i] = subs(pr = 0, f[i](pr)), i = 1 .. 2)};

Error, invalid input: f[1] uses a 2nd argument, pr, which is missing

Correct yourself the definition of Cs or f[i]'s before going further.

Have you ever considered that attaching an article might give the impression that you are asking the community to code and solve the equations it contains? Perhaps it is not surprising that no one then replies to you..

This put aside the error message is very clear: if a bvp problem is set over a compact subset of IR, let's say the interval [-1, +1], the the boundaries conditions have to be imposed at y=-1 and y=+1 and not aa some third point inside [-1, +1].
So your problem IS NOT a bvp problem.

This said there are essentially two ways to tackle it

1. Solve independently two bvp problems:
   1. One on interval [-1, 0] with unfixed values at y=0
   2. One on interval [0, 1] with still unfixed calues at y=0... but the same of the previous problem
   3.  Connect the left and right solution in order that your "internal conditions at y=0" are fullfilled.
   4. To do this you will have to solve some minimization ptoblem.
       
2. Solve a single bvp problem (the optionI choosed) over [-1, +1] with some bc unfixed and iterate over the solutions until your "internal conditions at y=0" are fullfilled.

The attached file explains you how to do this.

Last point: in OdeSys you have not given the parameter R some numeric value. So I arbitrarily choosed to take R=0.
If this choice does not suit yoy feel free to take another value.

The attached worksheet does not browse all the values of j as, for questions of time, I only considered j=1.
Here again it's up to you to complete the job if you agree with it.

symetry_paper_work_sand15.mw

see  plot/options

Q_figure.mw

Roughly speking you want to solve, with respect to 1 single indeterminate a__12,  33 polynomial equations in a__12 whose degrees range from 0 (!!!) to 3.
Are you seriously expecting to find a solution unless, maybe, the trivial one a__12 = 0?

Read this worksheet  help-parameter_sand15.mw

You will see that only 16 equations out of 29 (because 4 of the initial 33 equations simply do not contain a__12) have only in common the root a__12 = 0 and that the 13 remaining ones may have the same root provided the remaining parameters verify some relations.
If  by chance it happens to be the case, then the unique common root to all the 29 equations is a__12 = 0.

One you have understood what the attached file contains I advice you to reformulate what you want to achieve. Maybe fixing constraints on the remaining parameters or some other stuff.

The Warning, solutions may have been lost  comes from the fact that it is difficult for Maple to provide a general answer which depends on all the 7or 8 remaining parameters.

So the "solution" you expect is in fact an extremely complex piecewise structure, all the more complex to read that each clause and ppiece may have a lengthy expression.
In my opinion it is always a good idea in such a case to "synthetize" the expression you want to solve before really trying to solve it.
I suggest you this: Q_on_solve_sand15.mw

For two indeterminates named x and y:

restart
S := N -> ListTools:-Reverse([op(sort(add(seq(seq(x^i*y^j, j=0..N-i), i=0..N))))]):
R := S(2);
        [1, y, x, y^2 , x*y, x^2]

# EDITED:
assign( op([seq(p[i-1], i=1..nops(R))] =~ R) );


S(3); 

     [1, y, x, y^2 , x*y, x^2 , y^3 , x*y^2 , x^2*y, x^3 ]

For something more general you could use this  general_procedure.mw

Regional_plot_sand15.mw

Figure based on plots:-inequal:

Q_NEW_PLOT_sand15.mw

Plot 1 (for Plot 2 look inside the attached file)


Plot 1 (for plot 2 look inside the attached file)

Convolution is useful to build the PDF of the sum of two (or more) independent random variables but is of no help to compute the PDF of a function ( |.| for instance) of a random variable.

Furthermore, you write "..., but the "rules" of convolution, particularly setting the upper and lower boundaries of the integrals confuse me.".
This is very simple: as you are dealing with PDF of continuous random variables (so functions defined over R) the lower bound is minus infinity and the upper bound is plus infinity, that's all.
You can of course be more astute while integrating over the support of one of the two random variables, but I do not recommend you doing so: read Convolution_bounds.mw to see what mistakes you could do.
In fact, integrating over the smallest necessary range presents an interest only in numerical integration, but if you integrate exactly this trick has absolutely no utility.



Here is the solution in three steps.

Download DensityFunction_sand15.mw


BTW: in Statistics the distance you use is named Manhattan distance or "Taxicab distance".

 

Answers to questions 1 and 2:

# Question 1:
# Option gridlines = true/false manages the presence/absence of gridlines
# More details in plot/axis help page

# Question 2: compare this
z := x*cos(y):
p := plots:-contourplot(z, x=0..1, y=0..Pi, filled=true):
plots:-display(p);

# to this
plots:-display(
  plottools:-transform((x, y) -> [y, x])(p)
)

Question3: LOOK HERE or THERE for instance, and more generally use the search engine of this site.

My explanation intuition

On today machines 1 hyperthreaded proc corresponds to 1 phyisical proc + 1 logical proc.
In Grid language this means 1 core = 1 physical node + 1 logical node, so setting 4 nodes means 2 cores (=2 proc)

BTW: there is a course about Grid package somewhere here (a serie of two posts I wasn't capable to put the finger on,  but I'm sure I read them years ago...)
 

@brian bovril 

I've Googled the Round-Robin method (and discovered it was what is named the  "méthode du Tourniquet" in French, see picture below)

(this is my professional login, aka mmcdara at home)

Nevertheless I didn't find anything a Frenchie can easily understand.
So here is the last state of my work

RoundRobin.mw

Result:

L := [3, 4, x, x^2, x^3, sin(x)]:
select(t -> t::numeric and t > 1, degree~(L, x))
                             [2, 3]

# All the exponents but 0 and 1
L := [3, 4, x, x^2, x^3, sin(x), x^(-4)]:
select(t -> t::numeric and not(member(t,{0, 1})), degree~(L, x))
                          [2, 3, -4]

This piece of code "extracts" all integer exponents of x but 0 and 1

L := {3, 4, x, x^2, x^3, a(x^5), x^(-4), b(a(x^6)), exp(x), b(x^7, x^8)}:

# F := select(t -> hastype(t, `^` ), L); # is redundant
F := copy(L) # is enough because op discards 3, 4, x

while indets(F minus {x}) <> {} do
  F := op~(F);
end do:
F minus {x}
                   {-4, 2, 3, 4, 5, 6, 7, 8}

In case you would like to keep exponent 1 too:

L := {3, 4, x, x^2, x^3, a(x^5), x^(-4), b(a(x^6)), exp(x), b(x^7, x^8)}:

d := select(t -> t::numeric and t <> 0, degree~(L, x));
F := copy(L);
while indets(F minus {x}) <> {} do
  F := op~(F);
end do:
F minus {x} union d
               {-4, 1, 2, 3, 4, 5, 6, 7, 8}

Finaly, if some exponents are not integers and/or L contains terms of the form something*x^exponent you could use this code:

L := {3, 4, x, x^2+a*x^s, q*x^3, a(x^5), x^(-4), b(a(x^6)), exp(x^(-3/2)), b(r*x^7, x^p)}:

E := select(t -> t::numeric and t <> 0, degree~(L, x)):
F := select(t -> hastype(t, `^` ), L):
while F <> {x} do
  F, R := selectremove(has, F, x);
  F    := map(t -> if t::`*` then select(has, [op(t)], x)[] else op(t) end if, F);
  E    := E union R;
end do:
E;
                /                      -3      \ 
               { -4, 1, 2, 3, 5, 6, 7, --, p, s }
                \                      2       / 

Get_All_Exponents.mw

a := 3:

verify(a, 2..3, 'interval');                    # default means (2, 3)
verify(a, 2..3, 'interval'('closed'));          # here the interval is [2, 3]
verify(a, 2..3, 'interval'('closed'..'open'));  # interval is [2, 3)
verify(a, 2..3, 'interval'('open'..'closed'));  # interval is (2, 3]
                             false
                              true
                             false
                              true

You have to fix this big problem before

# Order of eq1 
D_terms  := select(has, indets(eq1, function), D):
Eq1Order := max(degree~(eval(eval(%, f = (eta -> exp(k*eta))), eta=0))); 
                               3

# Total number of IC/BC you impose
ficbc := select(has, {ics, bcs}, f):
NumberIcBc := numelems(%);
                               4

The number of IC/BC is not equal to the order of eq1.
Same thing happen to eq2: a 2nd order ode with 3 IC/BC.

Once you have fix these two problems you will be able (or maybe not...) to solve your differential system:

# Example 1 : eq1 with 3 ICs and eq2 with 2 ICs
# 
# As you then have "free" parameters [a, n, m] use option "parameters".

sol_1 := dsolve(
           {eq1, eq2} union select(has, {ics}, {f, theta})
           , numeric
           , parameters=[a, b, m]
         )
proc(x_rkf45)  ...  end;


# Example 2 : eq1 with 2 ICs + 1 BC and eq2 with 2 ICs.
#
# In this case you solve a coupled bv problem and you cannot use the option "parameters".
# You must instanciate [a, b, m] BEFORE calling dsolve, for instance:

data := [a=1, b=2, m=3]:

sys   := eval({eq1, eq2} union {f(0)=0, D(f)(0)=0, D(f)(10)=1, theta(0) = 1, (D(theta))(0) = b}, data):
sol_2 := dsolve(
           sys
           , numeric
         )

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

This error is quite common when you solve bv problems.
There exist several ways to get rid of this error but they are problem-dependent, in particular they may depend here on the values in data. So i will not waste tume to solve a problem wich could not be the one you are interested in.


Last point: writting things like Cp[f] := 4179 (for instance) is not safe (f is also the name of the function ). Look at this

restart
Cp[f] := 4179
                              4179
f := x^2;
                                2
                               x 
Cp[f]
                                 2 
                             Cp[x ]

Write Cp__f := 4179 instead (same thing for all the quantities whose names nontain [f])