Nice post!

What happens for r=r1 > 1680 in r2vsr1? Why this instability? What could be the interpretation for the value 1680?

@nm 

The problem is that sol is not a classic solution and Maple is not fully prepared to work with distributions such as Dirac.
(When using Dirac, the user is supposed to know what operations are allowed from a mathematical point of view).
In my opinion in such situations pdsolve returns a "formal" solution. It is the user who decides whether it suits his needs and in what sense the solution is to be interpreted (e.g. distributional). That's why pdetest cannot fully verify it.

BTW. The equation is very simple and the solution is simple too. Are you satisfied with it? Forget that it was found by Maple, do you accept it?

@Carl Love 

I was curious about the timings and I was a bit surprised that the seq version was about 10% faster.

@minhthien2016 

It was of course noted. As mentioned, all the listed solutions have the in-sphere centered at [0,0,0], So from the first solution ==> the tetrahedrons  {[-20+a, -3+b, -3+c], [-20+a, 3+b, 3+c], [20+a, -3+b, 3+c], [20+a, 3+b, -3+c]}
having the in-sphere at [a,b,c] in Z^3.

(Without the second condition the problem would be trivial.)

@ecterrab 

So, the help file says:

The eval option to subs also performs a full evaluation while keeping the semantics the same as for subs. This option should be used with care.

But your answer contains:
subs[eval] ... evaluate only the subexpressions where a substitution happened.

This explains the difference in my question. I think that the help file should contain this statement.
Now everything is clear. Thank you.

@ecterrab 

Thank you, the fix was really quick.
I still wonder about the difference between eval(subs(...)) and subs[eval](...) in my question.
Now, unlike value, the commands are builtin and unfortunalely the help page only mentions that the latter "should be used with care".

Happy New Year!

@Rouben Rostamian  

I think the behavior is intentional because it appears in an example (and in many help pages some essential information is found only in examples).
It is probably acceptable in the cases when the functions have always the same arguments. But othewise the results could be strange:

PDEtools[declare]();
diff(p(s,t),s) + diff(p(s,v),s) + diff(p(u,s),s);
     p[s] + p[s] + p[s]
lprint(%);
      diff(p(s,t),s) + diff(p(s,v),s) + diff(p(u,s),s);

@Christopher2222 

A default for an absolute width cannot exist because it should make visible the arrows for both vectors <0.5, 0.5>  and  <999, 999>.

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