Here is an idea how to achieve it over the integers (I think there was some
similar question lately)

Generate an strictly upper triangle matrix with integer entries. 

Generate an n-tuple from {+1,-1}. Take it, if it has an odd even number of "-1"s
and fill the diagonal. 

The determinant is +1. Replace exactly one arbitrary +-1 by +-2. Then det = 2.

Making it more arbitrary now generate elementary transformations invertible
over Z (like changing rows or columns [or transpositions]) and apply it.

I am not the one to talk about how "random" the result will be.

That is what I see opening you sheet. Very funny.

Int(Int(f(x)*g(y),x),y); expand(%);

Int(Int(f(x)*g(y),x=-infinity..infinity),y=-infinity..infinity); expand(%);

The following is from the cited book: read cos as the real part of exp and
find a path which reduces oscillation (I always admire such).

f:= x -> exp(ln(x)/x*I)/x;

'eval(f(z(t))*D(z)(t), z = 't -> t+I*t*(1-t)')';
simplify(%): # do not simplify it too much
g:=unapply(%, t);

Int(Re@g, 0 .. 1, method = _Dexp); evalf[40](%); # quickly found

              0.3233674316777787613993700879521704466510

That also works for crazy stuff like cos(ln(x^1)/x^2)/x^3, which may
serve as another test for MMA.

By simplification one can ignore k1, k2, but has R to be arbitrarly real.

a := cos(x)*exp((R*sin(x)-4*x)*I)*sin(x);

Plotting it for some data "shows" that the imaginary part is symmetric,

a = eval(a, x=-x); map(Im, %); simplify(%) assuming R::real, x::real;
is(%);

                                 true
The real part integrates to 0:

evalc(Re(a)); combine(%);
int(%,x=-Pi/2 .. Pi/2) assuming R::real;

                  cos(x) cos(R sin(x) - 4 x) sin(x)

         -1/4 sin(-6 x + R sin(x)) + 1/4 sin(-2 x + R sin(x))

                                  0

And the imaginary part can be found:

evalc(Im(a)); #combine(%); 
expand(%); #simplify(%);
2*Int(%,x=0 .. Pi/2) * I;

value(%) assuming R::real; #expand(%);
simplify(%, size);

It looks similar to MMA's automatic result, see at the end

For R = 2.3 I get -.693587157054838*I and for the integral

Int(a, x=-Pi/2..Pi/2);
subs(R=2.3, %); evalf(%);

                       0. - 0.693587157054844 I

with k1=1, k2 = 0 in the original formulation.

(4*R*Pi*(R^4-48*R^2+240)*BesselJ(1,R)+(36*Pi*R^4-480*Pi*R^2)*BesselJ(0,R)+(-2*R^5+224*R^3-1920*R)*cos(R)+(34*R^4-864*R^2+1920)*sin(R))/R^6

It is a bit ill-posed (made no attempt to correct it) and depends on Digits
as well. After pondering I passed to x = w^6 and normed the last variable.
The transformed function is named R.

Studying sign changes (and assuming continuity) there are lots of 'zeros'
for the suggested solution as well. I have some doubts in solutions using
double precision: "0" means a magnitude of 1E170 or so.

The last commands shows: changing the solution by the very last decimal by only 1 place
will cause the function to jump - in a magnitude of 10^171 around zero.

Download MP_some_Heun_zero.mw

I use cos(x^2-7) for a guess, the plot shows the idea: 

[f(x), sin(x^2-7)]; map(q -> diff(q,x), %);
plot(%, x=-6..6, color=[red, blue]);

I do it only for the positive axis (negatives are similar).

r:= k -> sqrt(7+(1/2+k)*Pi); # roots of cos(x^2-7) and a good guess

Then f1(r(k)) = -(-1)^k, hence the are infinitely many roots.

R:= k -> fsolve(0=f1(x), x=r(k)); # solver

seq(R(k), k=-2 .. 10); 

  1.49226378837861, 2.32016964600956, 2.92042014318168,
        3.41636824983574, 3.84873737108769, 4.23708504351411,
        4.59259995450624, 4.92238838366196, 5.23128723663552,
        5.52274710978730, 5.79929778739080, 6.06278493439976,
        6.31440696054959

If for a given x one wants the according k, x = r(k) then use

Z:= x -> round((x^2-7)/Pi - 1/2);  

Test it throug

[f(x), sin(x^2-7)]; map(q -> diff(q,x), %);
plot(%, x=100 .. 100.1, color=[red, blue]);

and 

xTst:=100.02;
0=f1(x); 
x0:=fsolve(%, x=xTst, fulldigits); #f1(%);
Z(%);
R(%);
f1(%);

It gives k=3182 and x0 = 100.025590094657

f1(x0) is not exactly 0, but that is due to limited exactness for
'rational' inputs, x0 is the 'best' approximation of a zero:

m,e:=SFloatMantissa(x0),SFloatExponent(x0);
'f1( SFloat(m-1,e) )': '%'=%;
'f1( SFloat(m,e) )': '%'=%;
'f1( SFloat(m+1,e) )': '%'=%;
                                                     -5
              f1(SFloat(m - 1, e)) = -0.5906311827 10

                                                   -5
                f1(SFloat(m, e)) = -0.2184335985 10

                                                     -5
               f1(SFloat(m + 1, e)) = 0.1537639854 10

Can not copy your 'pictures', but how about PDEtools[dchange] ?

Or eval(%, r = 's -> 1/u(s)'); ?

just post a new one - it is not a shame to have made errors - and for others it may be much more clear what is going on

I think you can reduce to show that

  8*EllipticPi((u-4)/u,1/(u+4)*((u-4)*(u+4))^(1/2))/((u-4)*(u+4))^(1/2) 

has Pi/2 as limit (that is the value of the integral)

"6/4"; parse(%);                                

            3/2

It is more Math than CAS ...

Any complex number z <> 0 can be written uniquely as r*exp(I*t), where 0 < r is real and -Pi < t <= Pi. So z is the very point on a circle in 0 with r = radius and t = angle, mathematical positive orientation (i.e. counter clockwise).

For t=0 you get the positive Reals, for t=Pi you get the negative Reals.

Multiplication of numbers is multiply the radius and add the angels (where you have to 'reduce modulo Pi').

For example squaring gives r*r * exp(I*t*2) [and reduce to the range]).

Now taking a 'n-th root' is just taking that root of the radius and devide the angle. N.B: for positive Reals one takes the positive root, which one needs to know for the radius.

For -8 = 8 * exp(I*Pi) hence you get 8^(1/3) * exp(I * Pi/3).

This is just rephrasing the help in terms of exp.

I think you can not do that in general, only numerically in some points.

You have 4 equations with 5 variables (and 2 parameters), and 4 of them
are exponents.

1. Generally having 1 variable more than equations defines a curve (or
something projected to the Reals which *may* be complicated).

2. Solving equations for exponents gives something periodic.

And in your case that additionally depends on 2 parameters ...


Writing p for phi[0] and with some manipulations your system is:

eqs := [2*((p-1)*q-p-1)*q*r^(-2+q)+(-2*n-2)*r^(-2+n)+(-3-a)*a*r^(a-2)+2*r^(-2+c), (-2+(2+(11+11*a)*a)*w^2)*r^(a-2)+(2+(4+2*n)*n*w^2)*r^(-2+n)-2*r^(-2+c)*w^2*(c+1), -2*(p+w+3/2)*(q+1)*q*r^(-2+q)-3*((p+1)*a+p-5/3)*a*r^(a-2), -2*q*(p-q-1)*r^q+(1+2*a)*a*r^a+(n+1)*n*r^n-r^c*c];
vars := {a, c, n, p, q};
params := {r, w}

The following gives what you want to see:

with(IntegrationTools);

Int(lambda(x)*(diff(y(x), x$2)), x);
Parts(%, lambda(x));

S1:-value(t = .2525,output=listprocedure);
f := rhs(op(3,%));
plot(f, 0 .. 1);
display(p2);

evalf(Int(f, 0 .. 1));
                          0.445408457258914

# check
numapprox:-chebyshev(f, x=0..1, 1e-8): 
#convert(%, rational):
eval(%, T = orthopoly[T]):
F:=unapply(%, x);

Int(F, 0 .. 1); value(%);
                            Int(F, 0 .. 1)
                          0.445408457329307

Likewise one can try to set up the DE for U = int( u(x,t), x )