Thanks to acer - he motivated me to give a another 'Maple answer'.

h := z -> (tan(abs(z)^(1/2))-abs(z)^(1/2))/abs(z)^(3/2), tan(x)=h(x^2)*x^3+x

x=Pi/4 means z = Pi^2/16 for z = x^2. Now use a Chebyshev-Pade approximation,
which is almost as good as a minimax approximation. Additionally I convert to
IEEE doubles by a function 'nearest' (which I once posted at the forum, one
can avoid it, if compiling the stuff). It is due to some experimenting to set
the digits for computing the approximation, 24 seems to be a good choice.

Digits:=24;                          # = 1.5* 16 decimals
chebpade(h(z), z=0..Pi^2/16, [6,8]): # new range (and instead of strange minmax)
evalindets(%, realcons, nearest):    # uses some non-intrinsic function, these
                                     # give special rational numbers = fractions
eval(%, T = orthopoly[T]):           # some stupid translation needed
convert(%, horner,z):                # may cause problems, see below 
evalf[18](%):                        # rationals to floats 
H:=unapply(%,z);
Digits:=18;

H := z -> 
 (.466433940928848490e-3+(1.11523043757015811+(-1.72989528300231501+(4.49729788979802074+
 (-.319393610864014770+(.545253820724829192e-2-.228858722931549036e-4*z)*z)*z)*z)*z)*z)/
 (.139930182278641728e-2+(3.34513159198135990+(-6.52796503942793992+(15.5613951451798833+
 (-6.34533642326600689+(.374356830875784296+(-.590832471659976894e-2+(.235408172559887276e-4+
 .461516094568396518e-9*z)*z)*z)*z)*z)*z)*z)*z)

approx:= x ->  H(x^2)*x^3+x;

That seems to have an absolute error of ~ DBL_EPSILON/2 and a relative error
of 1 ULP, without any warranty, but checked on sufficiently many points using
the above mentioned function 'nearest' - using IEEE doubles as inputs and as
outputs.

A better or even a proofed result would be quite challenging: note that already
the used Horner form for numerator or denominator is *not* correct w.r.t. to
numerical correctness in general (as far as I remember) and that problem will
increase with degree (so higher degree is not a recipe).

Also note, that one should look up and compare with Moshier's Cephes, SUN's
numerical libraries and Boot. For correctly rounded result there is crlibm
and others. Those will cover the reduction of the global Real function to
the interval as well.

Ignoring the 'b' that essentially asks for exponential function(x) = rational function(x).

Since the rational function has a unique solution (except its pole) I would expect that
this task has a periodic solution.

May be my feeling is wrong (sigh ...), but that 'tastes' like a classical question and
I would try to dig along that line (in ancient (?) literature/lectures).

Hm ...

Using continuous=true is nothing but using the anti-derivative, so it is a kind
of 'cheating' and only avoids doing that manually.

But actually Maple has the tools to find the discontinuities in that case, even
if it does not answer it for 1/20*ln(2*x^5+1):

Applying allvalues to RootOf would give it here - which is quite expensive and
would not work in general (beyond unit roots, 5 <= degree).

But "2*x^5+1; sturm(%,x,0,1);" would answer it quickly and restricting to the
case of interest - the intervall from 0 to 1.

May it would be worth to implement it for those cases, where the task is just a
polynomial problem over a real interval.

This may be complicated (and needs Maple installed, of course).

It is like calling Excel functions from C (without starting from Excel).

Why you want it and what are the details?

Generally it is not a good idea to inject Maple code (= 'create'),
but to have Maple code (already tested) to be called.

Example: you want to use Maple to provide a numerical value
(an integral or some special function). That will depend on the
settings (like Digits:= ...) and may results in spurious complex
values (due to numerical errors), while you expect real results
and only have these data type in C.

Since the numbers are small one can use floating point numbers. The following
ignores rounding issues for floats (by sqrt), but works well and using some
different printing it can be executed via 'evalhf' (N not too large):

check := proc(N::posint)
local A, a, b, k, d, e, big, small, g, cut, counter, cases, decide, minDist;
  counter := 0;
  cases := 0;
  cut := N*0.5;
  minDist := 3.0*N;
  print("a","b","k","d","e","g");
  A := Array(1 .. 5,datatype = float[8]);
  for a from 3 to N do
    for b from 23 to N do
      for k from b to N do
        for d from 4 to N do
          for e from 10 to N do
            counter := counter+1;
            A[1] := evalhf(sqrt(a^2+b^2));
            A[2] := evalhf(sqrt(k^2+d^2));
            A[3] := b+k;
            A[4] := e;
            A[5] := d+e-a;
            big := max(A);
            small := min(A);
            g := big-small;
            decide := g-cut;
            minDist := min(minDist,decide);
            if decide < 0. then
              cases := cases+1;
              print(a,b,k,d,e,g)
            end if
          end do
        end do
      end do
    end do
  end do;
  if 0 < counter then
    print("minDist =",minDist)
  end if;
  print("counter, cases = ",counter,cases);
  return NULL
end proc;

For N=30 it needs about 2 seconds for the ~ 5*10^6 combinations and shows, 
that minimum(big-small-15) = 1, i.e. big - small is at least 16. 

st:= time[real]():
  evalhf(check(30)): 
`time used` = time[real]() - st;

                     "a", "b", "k", "d", "e", "g"
                           "minDist =", 1.
                   "counter, cases = ", 571536., 0.
                          time used = 2.218

In 9.5 it is ok as well.

I can not believe that it depends on your version 14 (though bugs
are always possible), at least I first would look whether there is
an upgrade for 14.

Try to do it again with some "restart" before the integral.

Without the assumptions even in that old version I get a result
which can be read as "as far as A and B do not produce a pole".

PS: usually you simply can write "assuming 0 < A, 0 < B"

Too much printing may be a pain, especially in the standard interface, but anyway:

check:=proc(N)
for a from 3 to N do
  for b from 23 to N do
    for k from b to N do
      for d from 4 to N do
        for e from 10 to N do
A:=Array([sqrt(a^2+b^2), sqrt(k^2+d^2), b+k,e,d+e-a]); 
big:=max(A):
small:=min(A):
g:=big-small;
#print(whattype(g));
#printf("%g %g %g %g %g %g %g %g\n",a,b,k,d,e, g, big, small):
#  if (big-small)
if signum (big-small-15) < 0 then 
  printf("%g %g %g %g  %g %g \n",a,b,k,d,e, g);  
else
  end if:
end do:  #e
end do:  #d
end do:  #k
end do:  #b
end do:  #a
end proc;

Then the following terminates in ~ 10 minutes

st:= time[real]():
 check(30): 
time[real]() - st;

Persoanlly I appreciate that Maple stuff (never will like that word, those people are real persons, working) answers to questions.

But my attitude would be: " ... I forwarded it to Maplesoft Technical Support".

A user has a problem. Asks here. Gets an answer by an attentive employee actually saying "ask another employee".

That is not a good way. May be: internal note to Desk, which answers.

But that just my personal thought. Since it happens more than once it may be worth to ponder, a bit at least - as company.

Manually you can go this way (operating on functions, and evaluating them):

  Delta := h -> unapply(h(x+1)-h(x),x);

  Delta(f)(n);
                           f(n + 1) - f(n)

  (Delta@@2)(f)(n);
                     f(n + 2) - 2 f(n + 1) + f(n)

  (Delta@@20)(f)(n);


Likewise (if you know the very library)

  with(LREtools);

  delta(f(n),n);
                           f(n + 1) - f(n)

  delta(f(n),n,2);
                     f(n + 2) - 2 f(n + 1) + f(n)

  delta(f(n),n,20): sort(%); # somewhat longer

You may have should said that this is computing BesselK(n,z), n=0,1,2 and that
would also allow to judge results - they seem to be 'single precision'.

Found your stuff (as part of a Fortran program) on a Chinese server, which is
a bit 'strange', already its IP = 219.136.249.116 avoids to answer pings and
it offically is http://www.mesand.com providing pictures from flowers as the
only server on that IP.

Do you know who is the author (certainly older stuff)?

There are some issues and a bit can be said, even if you do not tell us about
the range of your x:=XIn (is it Real, Complex or an Integer? What magnitude
and what is the most likely range?).

Note, that for x=0 or x=5 it is not defined.

As Roman Pearce in the Maple newgroup already said you may use Arrays of certain
type and for hardware type it would give a lot. And would allow to compile, I guess.

If you do not use the global arrays outside your proc, then it may be faster to
have them inside (have not tested that)

Also avoid mixing datatypes: multiply by 0.5 instead of dividing by 2, dito XVal.
For C11 write ln(5) - ln(XAbs), division is expensive and ln(5) can be stored.

Then notice, that for 5 <= |x| your NST always equals 2. Thus you neither need to
compute it nor is there a need for 'add' in that branch.

Not sure, whether your C13 should be done at the end in that branch.

For the first branch |x| < 5 it certainly would speed up, if not using 'add'
4 times over the same indexing, I guess one loop would be faster (and would
allow to compile the code), the powers of XVal only have to be complied once
in that case and that can be done by xv_new = xv*xv_old instead of expensive
calls to the power function.

Also compute 1/x only once.

After that codegen[optimze](G, 'tryhadrd') should give you a good code I hope
(have not tried so far).

Finally your return is an Array of 3 entries, I think that is expensive.
I would say it is better to provide that as hardware float input argument
and fill in the results (you would use such for compiling anyway) and return
just nothing (or some status code).

If your x is very large then you may want to consider asymptotics.

1. Generally this should be put at the 'homework section', do not mind to classify at such

2. If you do not provide numerical data the you expect  an answer guessing them. Just
try, as gkokovidis suggests.

3. In case that this is a honst question: use "allvalues(%)" to get a general result and
provide your data after that, if needed.

4. If expecting trigs and exps: the SW does not know your desires, it has its own will.

my answer here is just a test ...

PS: could you let us know about it ?

Find attached my symbolic way, which gives your 1st plot, if being feed with data, so
it looks like a 'smoothed step function'.

I guess that Maple mixes up where Heaviside jumps, i.e. can not correctly find the zero
depending on assumptions and the integration bounds. 

One reason may be the inner integral, having the outer variable in its bounds. Changing
variables makes the inner integral to live on the unit interval (but the integrand is
more complicated by a linear transformation).

Then I evaluate the inner integral, but I prefer to use 'assuming' instead of 'assume',
see the sheet. Finally did the overall task (with some 'simplify/size' inbetween).

For feeding data I converted to rationals to avoid additional numerical problems.

Note that I replaced your ugly variable just by 'tau' and used another way to define
Heaviside(0).

MP_int_heaviside.mws

PS: so which of your plots is correct ?

The best would be, if you provide the special cases/values for which it fails.