Here's one way to do it. I'm going to name the procedure Apply because it applies its first argument to its second argument.

Apply:= proc(f, x)
   return f(x)
end proc;

g:= x-> x^2;
Apply(g, 2);

The above should return 4 of course. Make sure that it does that for you.

So, it's actually much simpler than what you were trying. Actually, the above is more elaborate than what's needed. It can be shortened to

Apply:= (f,x)-> f(x);

Let me take apart your code line-by-line to explain some of its problems. Your first line includes the parameter declaration var:= x. This makes var the parameter, and its default value will be x if SomeProc is called with less than two arguments. This is not a problem per se; I just want you to be aware of what your syntax means to Maple.

Your next line is f(var):= g. This is a very dangerous line, and I'm surprised and disappointed that the syntax checker allows it. It's dangerous because it allows the modification of something, f, which is external to the procedure and hasn't even been declared.

The command LargeExpressions:-Veil is pretty much tailor made for your problem. See ?LargeExpressions.

You mentioned using Newton's method, so I assume that you can do some iteration in Excel. I am about to present a very simple single-line procedure which when iterated exactly three times in hardware-float arithmetic will compute W(x) with a maximum absolute error of 6.7e-13  over your interval of interest, 0.3..30. The procedure is so simple that the only constants that it contains are 1, 2, and 3 (nary a decimal point!), and the only operations are 6 additions, 5 multiplications, 4 assignments, 3 divisions, and 1 natural logarithm. So that's only 57 operations over 3 iterations.

Here it is:

W:= (x,y)-> ((X,g)-> (S-> x*(1-(S+g/2)/(X*(S/g+1)+g/3)))(X^2))(x+1, x+ln(x/y));

"Whoa!" you say, "What do you mean by 'iterate' when the procedure takes more than one argument?"

It's technically called "folding"; I mean computing W(W(W(x,x),x),x). (I don't recommend that you evaluate that symbolically because the result is massive, as is typical with Newton's methods.) You can verify what I said about the absolute error with this plot:

plot(evalhf@((x-> W(W(W(x,x),x),x)) - LambertW), .3..30);

And this plot shows that the relative error is less than 2.7e-13 in hardware-float arithmetic:

plot(evalhf@((x-> W(W(W(x,x),x),x))/LambertW - 1), .3..30);

Now you say, "How did you derive that?" To which I say "with meticulous hand-craftsmanship, decades of practice, and a little help from Wikipedia and Maple." I'll give more details later. In the meantime, browse the Wikipedia article "Householder's method". My procedure W is a third-order Householder iteration. The fact that it's iterated three times is a coincidence.

I doubt that Excel will accept my terse functional-programming style, which isn't required for the efficiency. It's trivial to deconstruct, but I'll do it for you anyway:

W:= proc(x,y)
local g, X, S;
   g:= x+ln(x/y);
   X:= x+1;
   S:= X^2;
   return x*(1-(S+g/2)/(X*(S/g+1)+g/3))
end proc:

I agree that the infomation for each constant should be stored better. I think Record would be the best format. But even as it is, you do not need to know the position of derive in the sequence (which could vary anyway) in order to extract it. You can do

eval(derive, select(type, [ScientificConstants:-GetConstant(e)], `=`))

This will also act gracefully in cases where there is no derive.

AFAIK, circular derivations are not used in the package, and I'm glad for that.

Builtin procedures sometimes act weirdly when used in "naked" form, i.e., without (immediate) arguments. The most common case for me is eval. You can work around this by putting the builtin in a wrapper:

curry((e,i)-> `?[]`(e,i), f)([1]); 

Try this

funcCoeffList:= subsindets[flat](
   (Funcs:= [select(
      hastype, 
      indets[flat](
         subsindets[flat](expression, function, f-> One*f),
         `*`
      ), 
      function
   )[]]),
   {function, identical(One)}, 1
);
Funcs:= eval(Funcs/~funcCoeffList, One= 1);

In case it's not clear, I make One a dummy coefficient of all functions so that it'll be extracted properly, then I set its value to numeric 1.

Yes, Maple has a command nearly tailor-made for this job: Iterator:-TopologicalSorts (I don't know where that name comes from (and I've taken four graduate courses in topology)).  Here is the procedure code, and I'll provide some commentary after it. Warning, this is quite dense code, but it's super fast.

L:= proc(c::posint, L::And(list, satisfies(L-> nops(L)>=2*c)))
local p:= L[], A:= p[1+c..-c-1], P:= (p[1..c], p[-c..-1]); 
     [seq(
          [P[[seq(p[1..c])]], A, P[[seq(p[-c..-1])]]],
          p= Iterator:-TopologicalSorts(
               2*c, rhs~(op(3, rtable((1..c)$2, (i,j)-> i<c+j, 'storage'= 'triangular'['upper'])))
          )
     )]
end proc:

Examples of usage:

L(3, [a1,a2,a3, a,a, b1,b2,b3]);
       (output omitted)

nops(%); #Count entries of previous output.
       57

Many of your input parameters are redundant (which is dangerous in computer code). All that's needed is the value of c and the list of objects to permute that has at least 2*c entries. The list length m isn't needed. The as in the center are not needed, nor their count, and they don't even need to be the same. Also b is not needed. All that information comes from the input list itself.

Here's a speed test:

Perms:= CodeTools:-Usage( #command that measures time and memory usage
     L(6, [a||(1..6), b||(1..6)])
): #Very important! End this with a colon because the output is huge!
memory used=2.10GiB, alloc change=495.45MiB, cpu time=22.97s, real time=12.52s, gc time=13.78s
nops(Perms);
                            1566813

So, that's more than 1.5 million length-12 permutations in less than13 seconds.
 

Since it's such a common Question, it bears repeating that any number of plots with the same number of dimensions (2 or 3) can be combined with plots:-display as Kitonum has shown. However, for this common situation of showing a curve along with some related points, I prefer to do it with a single plot command, like this

X:= <0, 3.38, 5.21, 6.97, 8.4108, 10.099, 10.9232, 11.8091>:
Y:= <0, 0.760e-1, .1275, .1994, .2286, .3222, .3637, .999>:
h1:= solve(Vdc = 0.15e-2*sqrt(2.53669508e8*u^3-6.06101011e8*u^2+3.46343435e8*u), u):
nC:= nops([h1]): #number of branches
plot(
   [h1, <X|Y>], Vdc= 0..11.5, thickness= 0, color= [magenta$nC, blue], 
   style= [line$nC, point], symbol= asterisk, symbolsize= 24
);

The CDF of NegativeBinomial(r,p) can be expressed in terms of an analytic function, the regularized incomplete beta function I_p(r, x+1). This function isn't directly implemented in Maple, but it can be expressed as an integral "regularized" by dividing by an ordinary Beta function. So,

Statistics:-CDF(NegativeBinomial(r,p), x) = 
   Int(t^(r-1)*(1-t)^floor(x), t= 0..p)/Beta(r, floor(x)+1);
NB_CDF_1:= unapply(rhs(%), [r,p,x]):

A relevant identity is I_p(r, x+1) = 1 - I_1-p(x+1, r), which seems to work better for finding large quantiles for small p. So, we also have

NB_CDF_2:= (r,p,x)->
   1 - Int(t^floor(x)*(1-t)^(r-1), t= 0..1-p)/Beta(floor(x)+1, r);

You may confirm that either of these functions produce CDF plots identical to what you showed.

plot(x-> NB_CDF_1(1.965, 0.5, x), 0..10, view= [default, 0..1]);

The quantiles can be found by applying a continuous root finder to either of these CDFs. As I said, I found the second worked better in the region of interest. I used this as my quantile function:

NegativeBinomialQuantile:= proc(
   r::nonnegative, #"number of failures" parameter
   p::And(positive, satisfies(p-> p < 1)), #probability of success of each trial
   P::And(positive, satisfies(P-> P < 1))  #percentile
)
local t;
   Digits:= max(10, Digits);
   #Use common sense! Remove `floor` from strictly monotonic
   #functions passed to continuous root finders. Instead, use
   #`ceil` after root is found.   
   ceil(RootFinding:-NextZero(
      Q-> (1-P)*Beta(Q+1,r) - int(t^Q*(1-t)^(r-1), t= 0..1-p, numeric),
      0, guardDigits= max(5, 15-Digits), maxdistance= 10^Digits, _rest
   ))
end proc:

And I used it to make this plot:

plots:-loglogplot(
   p-> NegativeBinomialQuantile(1.965, p, 0.98), 1e-5..0.5,
   labels= [p, `98th %tile`], labeldirections= [horizontal,vertical],
   title= "The 98th %tile of NegativeBinomial(1.965, p) vs p",
   numpoints= 49
);

The best place to get help with Maple is right on this forum, MaplePrimes. So, ask a Question. I'd guess there are about a hundred plotting commands. Probably the best to start with is plot:

Example:
plot(x^2 + 1, x= -2..2);

The command is implicitdiff. Note below that I changed your assignment to w into an equation, and I added a multiplication sign to your second equation:

restart:
FuncRels:= { #Define implicit functions via set of equations:
   w = exp(exp(x*b)*(r-1)/(1+varphi*exp(x*b))), 
   ln(r) = varphi*exp(x*b)*(r-1)/(1+varphi*exp(x*b))-1
}:
Funcs:= {w,r}: #Define which *two* are considered dependent variables.
dwdb:= Diff(w(b,varphi,x), b) = #This line is just an inert label.
   implicitdiff(FuncRels, Funcs, w, b); #compute dw/db
dwdvarphi:= Diff(w(b,varphi,x), varphi) = #This line is just an inert label.
   implicitdiff(FuncRels, Funcs, w, varphi); #compute dw/dvarphi

The hard part is deciding which two are the dependent variables. Obviously w is one, but why is r the other? The only answer that I can give to that is that I tried {w,x} and I didn't like the results.

The lines that I called "inert labels" are strictly for your future (programmatic) reference. You can change them to whatever you want or eliminate them entirely.

A note on your terminology: This is called finding an (implicit) derivative. "Solving a differential equation" refers to solving an equation (or set of equations) that already has derivatives in it to obtain expressions for the dependent variables (or functions).

I was the author of the earlier Answer that you linked. Yes, that algorithm was optimized for and only possible for p=2.

Here's an implementation of Lucas's theorem for any prime p:

`mod/Binomial`:= proc(N::nonnegint, K::nonnegint, p::prime)
description "Computes binomial(N,K) mod p, p::prime, by Lucas's theorem";
option
   author= "Carl Love <carl.j.love@gmail.com> 30-Oct-2018",
   reference= "https://en.wikipedia.org/wiki/Lucas%27s_theorem"
; 
local r:= `if`(K > N, 0, 1), n:= N, k:= min(K, N-K);
   while k > 0 and r > 0 do 
      r:= r*irem(binomial(irem(n,p,'n'), irem(k,p,'k')), p) 
   od;
   r mod p
end proc:

Example of usage:

Binomial(100000, 50000) mod prevprime(50000);
      4

If you know that the number of boundary conditions (BCs) required for the system should be n (because you know how to count the differential order) and the error message says that you need more than n, then the cause is an extra symbolic parameter (symbolic = not numerically specified). This is often because you've spelled a parameter name different ways in your input.

You can list, in order^[*1], all the "top-level"^[*2] names used in any Maple expression expr (including a container expression such as a list or set of differential equations) by

indets(expr, name); 
lprint(%);

(The command lprint(%) simply makes the output of whatever immediately previous command easier to read "with a fine-tooth comb", and it's not necessary if the first command already makes what you're looking for obvious.)

The reason that dsolve complains the way that it does (rather than just saying "unspecified parameter found") is that you can have dsolve(..., numeric) solve for unspecified parameters of ordinary differential equation (ODE) boundary value problems (BVPs) by giving one extra BC for each such parameter.

Read my recent Post "Numerically solving BVPs that have many parameters". Just concentrate on the first few sections that cover the inputting of the BVP system. The later sections cover elaborate plotting techniques that may be "too much information" for you at this point.

^[*1]The order is fixed and chosen by Maple, but it usually bears some resemblence to alphabetical order.

^[*2]Names that only occur in indices are not "top-level" and are not listed seperately. Since they're listed with their "parent" names anyway, this is usually not an issue of concern. (Example: For t[final], final is the index and t[final] is the parent. The indets command that I gave would only list the latter.)

The bounds used to plot a function of two variables (by any of the numerous commands that do that) can usually be specified as the limits of integration of the corresponding iterated (double) integral, as in Kitonum's answer. This remains true in alternate coordinate systems: spherical, cylindrical, etc. I find this knowledge to be a great help to "solidifying" my mental visualization of 3-D problems.

It can be done like this:

subsop((indexList=~ ())[], L);

However, this is an inefficient operation for long lists, so if you find yourself needing to do this often, consider replacing lists with tables or Vectors.