Premature calculation occurs. Use  unapply  instead of  -> :

restart:
p := unapply(sum('a__||k' * x^k, k=0..5), x);
p(x);    # returns a0+ a1x +  ... + a5x5
p(1);   # returns a0+ ... + a5
p(0);   # returns a0

I removed the multiplication sign after  add  in  add*(d(k, j)*(diff(theta[j](t), t, t)), j = 1 .. 7) . Now everything works correctly:

restart;
for k from 1 to 7 do 
tau[k] := add(d(k, j)*(diff(theta[j](t), t, t)), j = 1 .. 7)+add(add(c[i, j, k]*diff(theta[i](t), t)*diff(theta[j](t), t), j = 1 .. 7), i = 1 .. 7)+phi[k]; 
end do:

Try these examples:

with(plots):
with(plottools):
with(CurveFitting):
BS:=BSpline(2, x);
plot(BS, x=-1..3);
op(BS);
op(6,BS);
plot(%, x=1..2, view=[-1..3,0..1]);
 

First click on the big letter T on the instrument panel:

Addition.  What is your goal? For a presentation, this is entirely appropriate, but if you want to perform various mathematical manipulations with this double equality, it is better to write it as a set or a list of two equalities. Then, for example, you can easily plot this line, check whether the point lies on this line, etc.

Examples:

L:={(x-2)/(-3) = y+4, y+4 = (z+5)/2}: # Specification of the line as a set of two equalities
P:=[x=2, y=-4, z=-5]: # Specification of a point
is(`and`(eval(L,P))[]); # Check that the point P lies on the line L
plots:-intersectplot(L[], x=-5..5,y=-5..5,z=-5..5, thickness=2); # The plot of the line L

Edit.

in Maple 2017.3:
 

Download simp1.mw

Edit.

solve command, which is part of the core of Maple, is one of the main and frequently used commands. It is unfortunate that it can not cope with such a simple equation. I think that developers should deal with this situation.
We can still solve this equation with solve command in two steps, first solving for  sin(t) :

solve(exp(2*sin(t))-1=0, sin(t));
solve({sin(t)=%, t>=0, t<=16}, allsolutions, explicit);
                                                            0
                    {t = 0}, {t = Pi}, {t = 2 Pi}, {t = 3 Pi}, {t = 4 Pi}, {t = 5 Pi}
 

If you do not know which command to use to finally simplify your numerical expression, it is useful to first calculate it approximately accurately (by  evalf  command), and then, for complete certainty, select the appropriate command for symbolic simplification, in fact already knowing the result.

S:=-1/2*((-30*sqrt(3)+81*I)^(2/3)+21+2*sqrt(3)*(-30*sqrt(3)+81*I)^(1/3)-6*(-30*sqrt(3)+81*I)^(1/3)+I*sqrt(3)*(-30*sqrt(3)+81*I)^(2/3)-21*I*sqrt(3))/(-30*sqrt(3)+81*I)^(1/3);
evalf[20](S);
factor(S);  # or
radnormal(S);
 

Your set is just one octant in 3d (in fact, all points in the space with negative coordinates) .

Download PS.mw

In this example, we generated an arbitrary polynomial of the third degree in 5 variables. But the method itself is of a general nature and it can be used to generate a polynomial of any degree and any number of variables.

polynom.mw

Edit.

Here is a simple procedure that solves the problem for not too large  d  and  K :

Pell:=proc(d::posint,K::posint)
local k, y, x, L;
k:=0; L:=table();
for y from 1 to K do
x:=isqrt(d*y^2+1);
if x^2-d*y^2=1 then k:=k+1; L[k]:=[x,y] fi;
od;
convert(L,list);
end proc:

Example of use:

Pell(27,1000000);

Output:
                  [[26, 5], [1351, 260], [70226, 13515], [3650401, 702520]]

Addition. Of course, much more efficient code can be written if we use  isolve  command, but I think that for a beginner the above will be enough.

Edit.

Specify the range for  x  variable explicitly:

with(DocumentTools):
A := [plot(sin(x)/x, x = -2*Pi .. 2*Pi, legend = sinc(x)), Sum(f(i), i = 1 .. N), plot3d(y^2*sin(x), x = -Pi .. Pi, y = -2 .. 2)]:
Tabulate(A, width = 60):

 

Download Tab.mw

Maximize (or Minimize) command only guarantees a local extremum, not a global one (see help for details). To get a global maximum, first plot a graph, and then narrow the range for the variable  x , that contains the maximum:

Addition. In order to immediately get the global extremum, you can use DirectSearch  package, which should be downloaded from the Maple Application Center and installed on your computer

https://www.maplesoft.com/applications/view.aspx?sid=101333

Download Max.mw

First you need to open your file, then  ctrl+ p

Another way: menu file -> print  or  print preview

Probably a double for loop is meant for summing such terms, if  i  runs from  0  to  m  and  j  runs from  0  to  n . A simple procedure named  Bxy  solves the problem:

Bxy:=proc(m::nonnegint,n::nonnegint)
local S, i, j;
S:=0;
for i from 0 to m do
for j from 0 to n do
S:=S+binomial(m,j)*binomial(n,i)*x^i*(1-x)^(m-i)*y^j*(1-y)^(n-j);
od; od;
S;
end proc:

Example of use:

Bxy(3,2);

                  

Here is the numerical solution:

restart;
pde:=diff(u(x,t), t$2)+diff(u(x,t), x$4)=0;
bc:=u(0,t)=-12*t^2, u(1,t)=1-12*t^2, D[1,1](u)(0,t)=0, D[1,1](u)(1,t)=12;
ic:=u(x,0)=x^4, D[2](u)(x,0)=0;
sol:=pdsolve(pde, {ic,bc}, u(x,t), numeric);
sol:-plot3d(t = 0 .. 1, x = 0 .. 1, axes = boxed);