If you want that eq and want to edit it: why do you not look for some Math editor / add in for MS Word ?

I changed your inital data to 'exact' numbers (and l1=r, l2=s to keep it readable
and put it into a reasonable sheet, not the strange standard *.wm).

After looking closer one can reduce the problem in 2 steps: first solve for x and y
with Eq 3 and Eq4.

Use that to feed Eq5 to express s=l1 by r=l2 and lambda.

Now you are at Eq1 and and Eq2, but those are 'only' in r and lambda.

These are 'algebraic' equations in r and lambda, r enters linear and 'solve' can
provide a solution - there are 6 of them (some are conjugated I think).

Find a sketchy sheet appended, hope you can finish your task from that.

MP_nonlinear_4variab.mws

solve for x or r?

degree in x it is = 16 and in r it is = 15

PS: whatever ... I doubt it can be done. Assuming it is x then divide out
((r*x^2+x-1)*(r^2*x^2-r*x-r+1) ) giving degree = 12. Now setting r = 1
Maple tells it is irreducible.
Otherwise set x=1 giving degree = 12 in r and seeting x=1 it turns out
to be a square of degree 6 and that is irreducible as well

I think it is an affine linear subspace in Z^4 of dimension 2 a linear subspace in Z^4 of dimension 3, find a sheet appended

MP_integer_subgrou.mws

Writing the assumed solution as formal power series 

A(x,y) = Sum( Sum(a(m,n)*x^m*y^n, n=0 .. infinity), m=0..infinity)

then I think one has to study 

(-x-y+1)*A(x,y) = quad + higher,
quad:=   a(0,0)+a(0,1)*y-a(0,0)*y+a(1,0)*x-a(0,0)*x-a(0,1)*y^2-a(1,0)*x^2
higher:= x*Sum(a(0,n)*y^n,n = 2 .. infinity)+
         y*Sum(a(m,0)*x^m,m = 2 .. infinity)

Note that quad has no mixed term x*y and higher also gives restrictions
for mixed degrees. I have not worked it out further. And you have to
give your initial 'values'. Perhaps even the motivation for the task.

I consider the two problems (original one and power series not starting in 0) as serious, since they are about hypergeometric summation (i.e. the quotient of coefficients is a rational polynomial in the summation index).

And I would expect that Maple can handle that correctly.

Has anybody already located the reason for that bug?

Sum(1/(6*k+1)!, k= n..infinity); 
                         infinity
                          -----
                           \          1
                            )     ----------
                           /      (6 k + 1)!
                          -----
                          k = n

value(%); evalf(%);

  1/6*3^(1/2)*exp(-1/2)*sin(1/2*3^(1/2))+1/6*3^(1/2)*exp(1/2)*sin(1/2*3^(1/2))+
  1/6*exp(1/2)*cos(1/2*3^(1/2))

                          0.673955774283463

Which is false: it does NOT depend on n, the starting index

You may want to provide all the values where m=0 or n=0.

Where I do not say, that Maple then would do the work (and 'genfunc' is a 1-dim function anyway)

Using an indicator function for the domain and Monte Carlo integration one should be
able to get some magnitude. Of course that's brute. And if the integration range is
too large the result will be too bad, here it is easy to see that |x| <= 1, |y| <= 1.

Int( (x,y,z) -> exp((x^2+y^2+z^2)^(3/2)) *
  # indicator for the domain, check for well defined expressions
  `if`(z <= (x^2+y^2)^(1/2) and 
       0 <= 1-x^2-y^2 and z <= (-x^2-y^2+1)^(1/2), 
    1, 0), 
  # range: 0 <= z, -1 < x < 1, -1 < y < 1  
  [-1..1, -1..1, 0 .. 1], 
  # do not expect to much from that brute approach, Monte Carlo may do it
  method = _MonteCarlo,
  digits=4); 
 
evalf(%);
                                2.545


Using the indicator as function to be called makes it slower

BTW the result does not become much better if using symmetries in domain and function
and avoiding expensive calls to 'sqrt'.

You may look up, where the Beta distribution is defined (Maple's help and Wikipedia?)
and whether your input is according to it: over the interval 0 .. 1. And the parameters
may have a direct interpretation, http://en.wikipedia.org/wiki/Beta_distribution, in
terms of your emprical data. Thus you can have a raw guess.

Similar for the (logarithmic) normal distribution the situation has to be checked first.


I just guess what you are doing: if these are pricing data Open-High-Low-Close then you
have to pass to (logarithmic) returns anyway, otherwise density estimation will make
no sense.

One can improve timings as follows:
 
1. Find a good initial guess for solving for the zeros. 
2. Find a representation for finding zeros which behaves well at the singularity z=0.


initial_x := z -> .166666*z+.9128709*(1.56332-2.*ln(z))^(1/2); 
g := (x, z) -> 3/5*x^2-1/5*x*z+1/60*z^2+
  ln(1/2*sqrt(2)*(1+erf(1/30*(6*x-z)*sqrt(15))))+ln(z*sqrt(Pi)) = 1/2*ln(30);;

s:= z -> fsolve(g(x,z), x=initial_x(z)); # to find the zeros


Then evalf(Int(z-> s(z)*PDF(Q,z), 0 .. 1)) returns .282975811112917 in ~ 0.5 seconds.

Find details attached for that kind of 'survival' problems.

MP_FsolveInt_surviva.mws
MP_FsolveInt_surviva.pdf

Additionally to Christopher:

If I have 32 bit programms (and there are many) then I install them into the folder
which is meant for them (something like "Program Files (x86)", it depends on your OS
language, MS is very idiotic and uses different names on command line mode and their
file explorer, if you have a non-english version).

Actually I use a different, new program folder without blanks (due to former frustration
that MS may not handle blanks correctly), like "programs_x86".

For that setting most old programs work.

But generally it is very hopefull that old programs work on OS for which they have
not been designed/compiled.

For example an old Maple 9.5 will not work correctly (and I gave up to look why).


Hope it helps, else you want to contact Maple's support via email.

As already said: 

I think one can not go the direct way, J = (N(t)*z - diff(N(t),t)*sqrt(30))/6 where
t = (6*x-z)/sqrt(30), N = cumulative normal distribution and thus it means to solve

'eval((N(t)*z - D(N)(t)*sqrt(30)), t = (6*x-z)/sqrt(30))';

               (N(t) z - D(N)(t) sqrt(30))|
                                          |    6 x - z
                                          |t = --------
                                          |    sqrt(30)

As z is in a small range around 0.5 that roughly means to solve Mills ration (or
its complement or similar?) to achieve a given value

N(t)*z / D(N)(t) = sqrt(30)

and I think one has to model that first (and not directly).

eval( N(t)*z / D(N)(t), t = (6*x-z)/sqrt(30)); lprint(%);

  (1/2+1/2*erf(1/60*(6*x-z)*30^(1/2)*2^(1/2)))*z*Pi^(1/2)/exp(-1/60*(6*x-z)^2)*2^(1/2)

eval( %, z=1e-3);                # Markiyan's test example
plot([%, sqrt(30)], x=-4 .. 4);  # showing ~ 3.5 as only solution

Edited: for that copy the output " (1/2+1/2*erf( ... " to plot it

As you already know r1,r2,r3 are roots as well as -3, 1, 4.

r2 < r3; is(%);
                                 true

r3 < r1; is(%);
                                 true

Thus it follows that r2 = -3 < r3 = 1 < r1 = 4


PS: I find it nice (and guess you implemented a solver using trigonometrics only).
Why dont you post your solver?

After some manipulation it can be written as hypergeometric 1F2

  1/4*1/a*hypergeom([1],[2, a+1],-1/4*z^2):
  % =convert(%, LommelS1); 

                                      2
                                     z
       hypergeom([1], [2, a + 1], - ----)
                                     4      LommelS1(a, -1 + a, z)
   1/4 ---------------------------------- = ----------------------
                       a                            (a + 1)
                                                   z

The limit for the LHS is 1/(4*a), a = 4.04 is your choice.

May be you also want to plot both for z = -1 ... 1