Logic:-Satisfy is new in Maple 2016. The variable p was forgotten from the local sequence.

As a workaround you can start with
local p
in your worksheet.

restart;
with(Logic):
local p;
p:=5;
Satisfy(x or y); p;
                     {x = true, y = false}
                               5


Edit. Now p is local to your program and the variable altered by Satisfy is :-p

To inspect the coordinates use plottools[getdata]

A:=plots:-implicitplot(x^2+y^2=1, x=-1..1, y=-1..1):
plottools:-getdata(A)[3];  # double click the summarized Matrix

convert(%,list,nested);    # if you really want a list

simplify(eval(x-y+4*z, [x = (u-v)^2, y =  u^2-3*v^2, z = (1/2)*v*(u-2*v)]));
      0

So, the surface is contained in the plane x-y+4*z=0.

(Of course x-y+4*z can be obtained via Groebner.)

r:=[(u-v)^2, u^2-3*v^2, (1/2)*v*(u-2*v)]:
with(LinearAlgebra):
E:=Vector(diff(r,u))^+ . Vector(diff(r,u)):
F:=Vector(diff(r,u))^+ . Vector(diff(r,v)):
G:=Vector(diff(r,v))^+ . Vector(diff(r,v)):

A := Matrix([[-diff(E,v,v)/2+diff(F,u,v)-diff(G,u,u)/2, diff(E,u)/2, diff(F,u)-diff(E,v)/2],
                     [diff(F,v)-diff(G,u)/2, E, F],
                     [diff(G,v)/2, F, G] ]):
B := Matrix([[0, diff(E,v)/2, diff(G,u)/2], [diff(E,v)/2, E, F], [diff(G,u)/2, F, G]]):
K := simplify((Determinant(A)-Determinant(B))/(E*G-F^2)^2);

                               0

You can use maxima. It is a free, complete and nice CAS which was ported to Android. It has many "common points" with Maple.
Note that on a phone or even on a tablet it is not comfortable to type longer code.

Mathematically this question is related to https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Arnold_representation_theorem

In your particular case, F(x,y) can be expressed as f(x/y) iff F is homogeneous of degree 0, which is not the case for F(x,y) = x+y-1.
Generally, F must have special properties in order to be expressed via functions of a single variable; e.g. if F has radial symmetry (i.e. F(u)=F(v) for ||u||=||v||) then F(u) can be expressed as f(r), where r = ||u||.

Your function F(n,N) = n+N-1 can be trivially expressed as z-1  for z=n+N, but obviously this does not help.

    floor((n+1)/2)

But why and how avoid the calculus of variations? Just pass the Euler-Lagrange equation and bc to dsolve.
Of course, sometimes we may use tricks.
E.g. for the provided example:

int( (f'(x) - 1)^2, x=0..1 ) = E - 1, so the minimum occurs for f' = 1 ==> f(x)=x.

If you have the ODEs

y'(x) = f(x, y(x))     (1)

y'(x) = g(x, y(x))     (2)

y'(x) = f(x, y(x))  + g(x, y(x))    (3)

then generally, the solution(s) of (3) cannot be expressed in terms of the solutions of (1) and (2), provided that such solutions exist.
Solving (3) using (1)+(2) is possible e.g. if f(x, . ) and g(x, . ) are linear. In this case (3) is also linear and its solutions can be expressed as y1 + y2, y1 being a solution for (1), y2 for (2).

It is difficult to debug 2d math code written by someone else. I don't use 2d math for input but I have noticed:

1. You used both p[n] and p__n. They are displayed the same. Please check. As you wrote it, your code needs both versions (not recommended!) but you must switch some of them (in u(x,t)).

2.  p__n ^ re should be Re(p__n). Same for im.

3. u(x,t) will depend on two variables. Use plot3d.

PS. It is not a good practice to use something like

a[n]:=n;
sum(a[n],n=1..10);

Here it will work. But add(a[n],n=1..10);  and   sum(a[k],k=1..10);  will fail.

Your integrals have a difficult removable singularity at 0 which is seen by some integration algorithms as a pole.

So, several methods must be tried and fine tuned.
I have set the method to _Dexp + minor changes and it works.

chebichef_cont-ok.mw

To inspect a large matrix without adjusting rtablesize and zooming, you can use the Browse Matrix tool.
Just double-click the placeholder matrix e.g.

A window with scroll bars will open and you will see all the elements.

You should use ApproximateInt only for didactic purposes (or fine tune its many optional parameters).

For example,

Student[Calculus1]:-ApproximateInt(sin(x),x=0..Pi);evalf(%);

    2.008248408

(the exact value being 2).

So, the reliable answers are given by evalf(Int(...))

Without epsilon this could remain unevaluated if the default accuracy cannot be obtained.
(this is because there is a removable singularity at 0 for your integrands).

To "prove" with Maple a nontrivial identity f = g, is is usually not enough.
In most cases, one defines h := f - g; and try to simplify to 0 using all sorts of acrobatics with simplify, convert, expand, assuming etc.

maximize(F,location);

works, as tomleslie stated.

Note that simplify(F) gives a better reprezentation.

Note also that your function is not continuous. If you use the almost equivalent definitition

G:=piecewise(t<0,0,t<=1,t^2,t<3,3-t);

then maximize(G,location) gives the same result

    2, {[{t = 1}, 2]}

but t=1 is not actually the maximum point, because the max is not attained.