@one man 

Yes, it works, but having direct access to the numerical procedures is usually easier and more flexible.

@Rouben Rostamian  

OK, I understand.  But at least if u0 is continuous, u is continuous everywhere (actually u(p) = ||p|| * u0(p / ||p||),  where ||.|| is the max norm (= infty norm), so it has a very simple form).

Edit. The problem has a solution. Using Perron's Method for the Dirichlet Problem, there exists even a harmonic u (in the interior of the cube) such that u is continuous at any continuity point of u0. Unfortunately a simple expression for u is not possible in this case; but we have uniqueness.

There are too many possibilities. There should be additional requirements.
For example, this simple construction works in any dimension. Let u0 be given on the faces.
For any p in the cube, p <> 0, let p' be the intersection of the half-line Op with the boundary of the cube.
Then an extension of u0 is given by u(p) := u0(p') * ||p|| / ||p'||,  u(0) := 0.
If u0 is continuous, u is continuous too. It's easy to give the analytic expression for u. 

@questerion 

So, your system has the form
Sum(Sum(a[i,j,k]*x[i]*x[j],i=1..n),j=1..n)=b[k];   # k=1..n

where a[i,j,k], b[k] are constants (k=1..n). If so, my remark applies.

You have no chance to solve your system this way unless most of the coefficients are zero.
For n unknowns x[i] and n quadratic equations the number of new unknowns s[i,j] = x[i]*x[j] is n*(n+1)/2.
Note that there are typically 2^n solutions for such a system.

@tomleslie 

Actually, det = 0. It appears to be small only due to roundoffs.

@Earl 

It's a minimal modification of Rouben's code.

restart;
with(plots):
r := <x(t),f(x(t))>:
v := diff(r,t):
v / sqrt(v^+ . v):
T := simplify(%):
N := < -T[2], T[1] >:
F := <0,-m*g> + n*N - mu*n*T*signum(diff(x(t),t)):
m*diff(v,t) =~ F:
DE := convert(%, list):
solve(DE, {diff(x(t),t,t), n}):
select(has, %, diff(x(t),t,t)):
de := %[]: 
myde := eval(de, {f = (x -> x^2), m=1, g=1, mu=1/10}):
ic := x(0)=0, D(x)(0)=1:
dsol1 := dsolve({myde, ic}, numeric, output=listprocedure):
pt:=eval(t, dsol1): px:=eval(x(t), dsol1):
tmax:=13.55:
parab:=plot([px(t),px(t)^2,t=0..tmax],color=blue):
animate(plot, [[[px(t),px(t)^2]],style=point,symbolsize=30], t=0..tmax, frames=100, background=parab);

@acer 

In the floating point case the difficulty is inherent because the rank is not continuous wrt the matrix. The Rank function should return FAIL for a matrix for which the mathematical rank is not continuous. For example, Rank should return 2 for the matrix <1.,2.;3.,5.> and FAIL (or maybe 1..2)  for <1.,2.;3.,6.> and should depend on precision. All these in an ideal world.

@Mariusz Iwaniuk 

You should care because:

J:=Int(convert((cos(2*x))/(1+2*sin(3*x)^2), exp), x = 0 .. Pi/2):
evalf(J) = evalf(value(J));

    -3.452793607*10^(-14) = .3318985824+.2373064033*I

@Axel Vogt 

@Rouben Rostamian  

Yes, you are right. The changes of variables failed for the intervals not included in [0,Pi].
It is sad that we have to check each simple step that Maple does :-(

P.S. Actually I was surprised that Preben's split seemed to work because I have tried myself similar methods and failed.

Edit.
We should choose only changes of variables with "unique" inverses (mainly for trig functions), Maple not being smart enough to find the correct branch.
Change(Int(f(x), x=3*Pi/2 .. 2*Pi), cos(x)=t);
         
is wrong. We'll have to use:

Change(Int(f(x), x=3*Pi/2 .. 2*Pi), x=2*Pi-arccos(t));

@Preben Alsholm

Yes, it's correct, simple, and can be done by hand (the 4 integrals cancel in pairs). Thank you.

My workaround was

T:=IntegrationTools:-Change(J, tan(x/2)=t):
subsindets(T, `=`, u -> (t=-infinity..infinity))/2:
simplify(expand(%)):
value(%);

I did not like it because the parity and periodicity is used and the degree of the polynomial is 12.
(It can be reduced to degree 6 with another change of variables but then evala() is needed.)

@tomleslie 

It is not correct because the antiderivative obtained by int(...,x) is not continuous.

@Rohith 

Yes. But this does not imply that such a function is unbounded; e.g. f : R --> R, f(x) = arctan(x).