@torabi 

It seems that you want to obtain some polynomial approximations for an ODE system. You impose e.g. K(0)=0 but in your system K appears at the denominator; it cannot work.

If you change
s2:=K(eta)(0) = 0.001: s3:=Omega(eta)(0) = 0.001:
then fsolve finds a solution, but I don't know if it is related to the solution of the ODE (provided it exists, which is problematic).

@Markiyan Hirnyk 

I don't think that you are in a good position to talk about seriouseness and also politeness here.

@Markiyan Hirnyk 

The post was edited. The Gamma stuff was in the second question of the OP.

@Markiyan Hirnyk 

This is not a scientific paper. The OP assumed that the sample comes from a Gamma distribution and I have shown how the parameters a,b can be obtained. That.s all.

plot/discont seems to be confused when f(x0) is undefined.
This is mathematically acceptable because the continuity makes no sense at a point which is not in the domain of definition.
Actually plot/discont seems to ignore the value of f at x0. Instead, it considers the values near x0, in particular the side limits (if they exist).
For example, f := frac ploted in [0,5] with option discont will have a vertical line at x=5 even if we define f(5):=1 (and so f is now continuous at 5 in [0,5]).

Note that
f:=x->piecewise(x<=4, frac(x), 1);
plots without a vertical line at x=4, ignoring that f(4)=0.

@jamunoz 

Of course the precision depends on the x-step but I don't know about an error estimate in computing the derivative [we would need for example  bounds for f''' ].
We can use higher-order methods for f'(x), e.g.

f'(x)  ~= (f(x-2h) - 4f(x-h) + 3f(x)) / (2h)

@eclark 

It gives more than that, namely the open ends of the interval; in this way tou can replace Open(b) e.g. by b - eps.
You can write a very simple proc to extract a,b and obtain the desired form.

torange:=proc(ineqs::{set,list}, var::name)
  local a,b,x, rr:=subs(var=x,ineqs);
  assume(op(rr));
  rr:=op([1,2],subs[eval]([Open = (a->a)], getassumptions(x)));
  eval(rr, RealRange = ((a,b)-> a..b));
end:

torange({theta <= Pi-arccos(-3/4+(1/4)*sqrt(13+16*sqrt(2))), arccos(3/4) < theta}, theta);

Expand what? The question is more than ambiguous (and your reply too, BTW).

@Matt C Anderson 

@Carl Love 
No, the modules are just fine. As I said, it was pure curiosity, but I don't use them very often (except packages) and I was trying to understand if in such situations they are really more efficient.
Using modules the debugging could be more difficult if something goes wrong.

Edit. You say: "the syntax of the inner procedure is checked when it's originally entered. But the assignment := is not done".
AFAIK, when defined, the body of the procedure is translated into internal code (that is what I understand by parse here). So, why should the inner procedure be "re-parsed"?
 

@Carl Love 

But in the procedure version of MakeTreeTable, when MakeTreeTable is defined (parsed), this implies (in my opinion) that all its body is parsed, including MakeTable.
Probably in other circumstances the use of modules is preferable.

@Carl Love 

Just curious: is there any benefit in using modules e.g. in MakeTreeTable versus the equivalent procedure?

MakeTreeTable:= proc(depth::nonnegint, width::nonnegint)
local
   i, node:=0,
   nextnode:= proc() node:= node+1 end proc,      
   MakeTable:= proc(depth :: nonnegint, width :: nonnegint)
      if depth = 0 then nextnode()
      else table([seq(nextnode()=thisproc(depth-1,width), i=1..width)]);
      end if
   end proc;
MakeTable(depth, width)
end proc;

@asa12 
Let us try to formulate mathematically the problem (if I understand it correctly).

You have a function space  FS, for example C^infinity(J)  (= the algebra of infinitely differentiable functions on the open interval J) . You also have an operator A defined (partially) on FS, in your case it is:

   A(a) = a' / (a * a' - a).

You want to find another operator f defined on FS such that

   f( A(a) + A(b) )  = f(A(a)) * f(A(b)),  for each a,b in FS.

Such an f is for example f = exp, or more generally, f(a) = exp(s*a),  s being a scalar.

Now you want an operator f, similar to diff. I suppose that this means:

  f( A(a) * A(b) ) = f(A(a))*A(b) + A(a)*f(A(b)).      [*]

(I suppose you want to exclude the trivial case f=0).

Maple probably cannot do this. And I am not sure if such an operator - other than f = D - exists. It may depend on FS, which you should know, because it is your ploblem. And I am not sure if [*] is what you want.

@asa12 

Maybe you want an Infinitesimal generator, see https://en.wikipedia.org/wiki/C0-semigroup#Infinitesimal_generator
But you must learn no ask the question in proper terms (preferably using some maths and without non-working Maple code), otherwise it is "quasi-impossible" to be understood.
 

@Carl Love 

Nice solution, vote up.
Don't you thing that using multi-indices is more efficient?
I wonder why the OP has no reaction.

@Kitonum 

It's OK but you should add the step about the equivalence, for those thinking that it is not obvious.

[ 2y=x(1-y^2)  <==> 2y/(1-y^2) = x,  because y cannot be +-1]

P.S. In my opinion, using Maple without knowing maths is not recommended.