Here is a procedure for this:

restart;
NestedSeq:=proc(Expr::uneval, L::list)
local S;
eval(subs(S=seq, foldl(S, Expr, op(L))));
end proc:

As an example, let's consider the numerological game "Cчастливый билет", known in Russia. A ticket in public transport is called lucky (a ticket is an ordered set of 6 digits, each digit is from 0 to 9), if the sum of the first three digits is equal to the sum of the last three. Below we find the list of all the lucky tickets:

LuckyTickets:=[NestedSeq(`if`(i+j+k=l+m+n,cat(i,j,k,l,m,n),NULL), [i,j,k,l,m,n]=~0..9)]: # List of all the lucky tickets
nops(LuckyTickets); # Total number of lucky tickets
seq(LuckyTickets[i],i=2..%,5000); # Examples of lucky tickets

                                                             55252
100100, 107701, 908791, 728782, 166373, 528654, 593935, 316226, 942807, 748397, 568388, 977689
 

This is impossible even if you know the function, because at the same value of the function, but at different points, the derivatives may differ. For example, cos(Pi/2)=cos(3*Pi/2)=0 , but the derivatives at these points are different. To calculate the derivative at individual points, it is better to use the differential operator  D .

In your example you can just write

diff(f(x),x)*eval(cos(f(x)), f(x)=0);

If you still want to apply  eval  to the whole expression, then the use of  D  operator is the simplest way:

eval(D(f)(x)*cos(f(x)), f(x)=0);  # Or
convert(eval(D(f)(x)*cos(f(x)), f(x)=0), diff);

Download Help_new1.mw

This can be done in several ways. Here are two:

with(plots): with(plottools):
A:=[1, 2, 3]: B:=[4, 5, 6]:
spacecurve(A+t*~B, t=0..1, color=red, thickness=2);
display(line(A, B, color=red, thickness=2));
 

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