@Carl Love 

Yes, I know, but the procedure is recursive and I can't do that unless the procedure is wrapped in a calling one which will have to increase Digits. I wanted to keep Cf very simple, otherwise some other things could be also added.
Notice also that Digits must be increased also for x (before entering Cf(x)).
After all, it is only a toy which is not going to be often used, don't you think so?

@Axel Vogt 

Yes, and I have mentioned this possibility. This happens iff x is in C \ (D union {0,1}) where C is the Cantor set and D is the set of the endpoints of the removed intervals i.e. D = {1/3, 2/3, 1/9, 2/9, 7/9, 8/9, ... }.

This set is "small", its Lebesgue measure being 0 (it is uncountable though!).
1/4 and 1/12 are in this set and their image cannot be computed exactly, even if they are rational.
I have also mentioned a:=Sum(2/3^(n^2),n=1..infinity) for which Cf(a) is irrational.

Edit. I will have to change Cf to its (initial in fact) version

piecewise(x<1/3, Cf(3*x,n-1)/2, x<=2/3, 1/2, 1/2+Cf(3*x-2,n-1)/2)

because in its actual state the approximation is < Digits for these cases. 

@Preben Alsholm 

But I am afraid that trying this for sin(x)=-1/2 the OP will be in a total darkness.

@AmirHosein Sadeghimanesh 

H:=(f,x,j) -> quo(f,x^(degree(f,x)-j),x);

@Axel Vogt 

You are right, thank you. BTW, Cf(1/2+eps) = Cf(1/2-eps) = Cf(1/3+eps) = 1/2 obviously :-)

@Axel Vogt 

The test value for x was actually x = 1/3 - eps/3.

As I mentioned, Cf should compute the exact value of Cf(x) using Cf(x,n). On my computer Maple failed (disconnected) probably due to insufficient memory.

But, Cf(1/3 - 1/(10^40000*Pi)) can be manually computed (using Maple, of course) exactly: it is 1/2 - 3/2^83839.

(Edited typo).

@Markiyan Hirnyk 

As I have explained, if you want 100000 digits, you must set Digits:=160000, such that Cf uses this precision.

If you want the exact value of Cf(1/3-1/(10^40000*Pi)) then call

Cf(1/3-1/(10^40000*Pi), n);

and increase n until the result is rational. But such huge values could exhaust the memory. (There is also the improbable possibility that the exact value is irrational; in this case it cannot be computed exactly with Maple).

If you want to test the precision I recommend to use values for x for which you know the exact answer (as I did).

My assertion concerning the exactness of Cf(x)  when it is rational is under the assumption that x is rational or symbolic such as 1/Pi  (and of course we count on Maple to compare correctly the irrationals by increasing Digits if necessary).
Of course, only rational numbers can be represented exactly by floats.

Now, if x is a float, then the precision decreases due to roundoff  errors.
After some experiments it seems that in order to have d correct digits for Cf(x) when x is a float, we should use Digits = floor(d*1.6).

Let's take an example.
a:=Sum(2/3^(n^2),n=1..infinity);
b:=Sum(1/2^(n^2),n=1..infinity);
# Then Cf(a)=b. [Note that a and b are irrational]

Digits:=1000*16/10;
fa:=evalf(a):
fb:=evalf(b):
Cf(fa):
evalf(%)-fb;   # ==>  -2.16*10^(-1013)  So, Cf(fa) is accurate to >1000 digits

Digits:=1000;
fa:=evalf(a):
Digits:=1000*16/10;
fb:=evalf(b):
Cf(fa):
evalf(%)-fb;   # ==> -1.04*10^(-637)

As you see, it is not enough to compute Cf(x) with floor(d*1.6). digits; the input float value x must also have the same precision!

@acer 

Maybe a good idea would be to have both a full ConditionNumber (according to the definition) and a (probably faster) ConditionNumberEstimate.

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))]);