Indeed as Carl noticed, int has problems because for the Cantor's function the set of points of non-differentiability is infinite. BTW the integrals int(Cf^n, 0..1) are known, and are expressed via Bernoulli numbers.

Cf is computed with Digits precision. Actually, for "most" x in [0,1], Cf(x) will be exact. If the result Cf(x) is rational (rather than float), then Cf(x) is "guaranteed" to be exact. Note that Cf(x) is rational a.e. (that is except in a set of Lebesgue measure 0).

@Markiyan Hirnyk 

I will come back later, I am busy now.

The Cf(1/P1) problem can be solved replacing type(x,numeric) with type(x,realcons).

@Markiyan Hirnyk 

For such a simple matrix having the condition number 6, I would accept say 5.9 for an estimation, but not 3 (!).
Unfortunatele it is very time-consuming for a user to localize the problem (or even impossible if compiled programs are involved).

@Markiyan Hirnyk 

It seems to me that the FPStruct is not correct (probably got confused by the unevaluation 'S' , ?).
v(n) = 0 satisfies the recurrence, but -x * sin(x)  + _C[0]*x^2  cannot be a solution.

@shani2775 

Why not this approach?

restart;
N:=8;
f:=unapply(add(d[k]*(x-a)^k/k!,k=1..N),x);    # d[k] = (D@@k)(f)(a)
g:=unapply(x-f(x)/D(f)(x),x);   # Newton's approximation
series(g(x),x=a); 
series(f(g(x)),x=a);

@mkonto 

The vector allSolutions depends on the variables t[1], t[2], ...
You must substutute each of them.
eval(allSolutions,[t[1]=7,t[2]=5]);

eval(allSolutions, t=0) works because 0[x] evaluates to 0; for t=1 it does not work because  1[x] remains unevaluated.

If you want a random vector from allSolutions, use e.g.

eval(allSolutions,[seq(u='rand()', u=indets(allSolutions))]);

@Carl Love 

In wikipedia b=0. But after a translation by a vector u_0 (where Au_0 = b), one obtains Rouben's formula.

@mkonto 

I have changed the procedure to work with Gradient without changing the Gradient vector to free.

with(VectorCalculus);
SetCoordinates(cartesian[x,y,z]):
n1 := Gradient(a*x^2+y^2+z^2):
nf:=ConvertVector(n1,free):
vf:=Proj(<<1|2|3>,<2|4|6>>, <1,2>, nf);
v1:=Proj(<<1|2|3>,<2|4|6>>, <1,2>, n1);


V1 and vf will be displayed as

But if Proj is defined after loading VectorCalculus, they will will be displayed:

@Markiyan Hirnyk 

Thank you.

@Dave L 

I did not noticed there the interdiction to post (and discuss) the result of a Maple command such as

showstat(int)

which BTW does not contain any copyright mention. After all, Maplesoft offered this possibilty, it is not a "reverese engineering" made by a user. Otherwise it would be problematic to post the result of int(x^2,x).

@MariaL 

Just inspect the Cayley table. Foget about 90, it was probably a typo.

@krismalo 

fsolve without location finds only U__C <0 and the problem needs >0. That is why {U__C=0..1000, X=0..360} was added, but fsolve failed (not unusual!).

_Z1 and _B1 are arbitrary integer and binary constants contained in the general solution (see ?solve,details) .  They must be specified if a numeric solution is wanted.

@zia9206314 

p is your polynomial. I converted it in 1D math. In the original 2D worksheet the name Nu was somehow distorted and I replaced it by N [but you may revert to Nu if you want]:

p:=-(1/416179814400)*(N^4*k1^4*k2^7+7*N^4*k1^4*k2^6+21*N^4*k1^4*k2^5+35*N^4*k1^4*k2^4+35*N^4*k1^4*k2^3+364*N^3*k1^3*k2^5+21*N^4*k1^4*k2^2+1820*N^3*k1^3*k2^4+7*N^4*k1^4*k2+3640*N^3*k1^3*k2^3-23744*N^2*k1^2*k2^5+N^4*k1^4+3640*N^3*k1^3*k2^2-94976*N^2*k1^2*k2^4+1820*N^3*k1^3*k2-118160*N^2*k1^2*k2^3+364*N^3*k1^3-22064*N^2*k1^2*k2^2+49168*N^2*k1^2*k2-1435648*N*k1*k2^3+24304*N^2*k1^2-2871296*N*k1*k2^2-1216192*N*k1*k2+9625600*k2^3+219456*N*k1+9625600*k2^2-3990528*k2)*N^3*k1^3+(1/147456)*(N^2*k1^2*k2^4+4*N^2*k1^2*k2^3+6*N^2*k1^2*k2^2+4*N^2*k1^2*k2+N^2*k1^2+56*N*k1*k2^2+112*N*k1*k2+56*N*k1-640*k2^2-640*k2+144)*N^2*k1^2-(1/14745600)*(N^3*k1^3*k2^5+5*N^3*k1^3*k2^4+10*N^3*k1^3*k2^3+10*N^3*k1^3*k2^2+5*N^3*k1^3*k2+120*N^2*k1^2*k2^3+N^3*k1^3+360*N^2*k1^2*k2^2+360*N^2*k1^2*k2-2944*N*k1*k2^3+120*N^2*k1^2-5888*N*k1*k2^2-1520*N*k1*k2+1424*N*k1-13312*k2)*N^2*k1^2+(1/2123366400)*(N^4*k1^4*k2^6+6*N^4*k1^4*k2^5+15*N^4*k1^4*k2^4+20*N^4*k1^4*k2^3+15*N^4*k1^4*k2^2+220*N^3*k1^3*k2^4+6*N^4*k1^4*k2+880*N^3*k1^3*k2^3+N^4*k1^4+1320*N^3*k1^3*k2^2-9344*N^2*k1^2*k2^4+880*N^3*k1^3*k2-28032*N^2*k1^2*k2^3+220*N^3*k1^3-21008*N^2*k1^2*k2^2+4704*N^2*k1^2*k2+7024*N^2*k1^2-205312*N*k1*k2^2-205312*N*k1*k2+14400*N*k1+409600*k2^2)*N^2*k1^2-(1/4)*(k2+1)*N*k1+(1/64)*(N*k1*k2^2+2*N*k1*k2+N*k1+4)*N*k1-(1/2304)*(N^2*k1^2*k2^3+3*N^2*k1^2*k2^2+3*N^2*k1^2*k2+N^2*k1^2+20*N*k1*k2+20*N*k1-64*k2)*N*k1+1;

Then the command works.
(It is much better to use 1D math because you can see exactly the content)

collect-expand-Nu.mw

@Markiyan Hirnyk 

I inserted some comments. Note that the code is short and does not use tricks, so I think that more comments are not cecessary.

@dorna01 

I have spotted an error in M2[2,1]; you have omega(1-25000000000000/...)
instead of (probably)  omega*(1-25000000000000/...)
but this does not change essentially the problem (i.e. the equation is almost 0=0).