Should be Pi instead of pi :

u:=(x,t)->Sum(sin(r*Pi*x/20)*(4/(r^2*Pi^2))*sin(r*Pi/2)*cos(r*Pi*t/20), r=1..1000);

plot3d(u(x,t), x=0..10, t=0..1);

We can see a certain periodicity, if we increase  t :

plot3d(u(x,t), x=0..10, t=0..100, numpoints=10000);

RealDomain:-eval(log[1/3](x)-log[sqrt(3)](x^2)+log[x](9), x=3^a);

Should be  add  instead of  sum  (or before eqs execute  i:='i': ). See corrected file   1_new.mw

Try

lambda:=unapply(5*Pi*sqrt(m^2/16+n^2/4),m,n): 
u:=(x,y,t)->0.426050*add(add(1/(m^3*n^3)*cos(lambda(m,n)*t)*sin(m*Pi*x/4)*sin(n*Pi*y/2), m=1..100,2), n=1..100,2):
plot3d(eval(u(x,y,t),t=0.01), x=0..2, y=0..2, style=patchcontour,axes=BOXED);
#  Or even easier
plot3d(u(x,y,0.01), x=0..2, y=0..2, style=patchcontour,axes=BOXED);

I replaced  infinity  with 100. Note also that I changed  sum  to  add . For finite sums (not symbolic), add  is preferable.

Edit.

If you have a polynomial  a*s^4+b*s^3+c*s^2+d*s+e  with symbolic coefficients  a, b, c, d, e  then you can factor it  (s-s1)*(s-s2)*(s-s3)*(s-s4)  explicitly (first finding its roots explicitly), but only the expressions for  s1, s2, s3, s4  will be generally very cumbersome. See example:

L:=[solve(a*s^4+b*s^3+c*s^2+d*s+e, s, explicit)]:
mul(s-p, p=L);   # Factorization in symbolic form
eval(%, [a=1, b=2, c=-1, d=-4, e=-3]);  # Replace symbols with numbers (in symbolic form)
evalf(%);  # Numerical factorization

See the example in  the corrected file  question_new.mw

See help on  fnormal  command.

Points:=[[2,4], [3,5], [4,7], [5,8], [6,11]];
Line:=plot(Points, color=blue):
points:=plot(Points, style=point, symbol=solidcircle, color=red, symbolsize=15):
Curve:=CurveFitting:-PolynomialInterpolation([[2,4], [3,5], [4,7]], x);
P:=plot(Curve, x=1..6, color=green):
plots:-display(Line, points, P, view=[0..7, 0..12]);

Addition. I did not understand the meaning of your question "And the codomain is also the set of positive integers?"

The points marked with a red cross are found as the points of intersection of your curve with the y axis, and the coordinates of the highest point by  maximize  command:

F:=convert((1.25*y-sqrt(abs(x)))^2 + x^2-1, fraction);   # Converting to a form with rational coefficients (for symbolic results)
solve(eval(F, x=0), y);    #  Red cross points
f:=solve(F, y);   # The top and bottom graphs  y  vs  x  (explicit equations)
maximize(f[1], x=0..1, location);    # The right green circle (exact result)
evalf(%);    # The right green circle (approximate result)

We can use the mathematical induction to get the general symbolic result, because

2^(n-1)+2^(n-1)=2^n

The symbolic answer is

f(n) = piecewise(n=0, 1, 2^(n-1))

Maple fails with this:

rsolve({f(n+1)=sum(f(i), i=0..n), f(0)=1}, f(n));

Another way is to use  PolyhedralSets  package:

with(PolyhedralSets):
ps := PolyhedralSet(convert(subs(sqrt(2)=evalf[4](sqrt(2)),[2*sqrt(2)*x-2*z+sqrt(2) >= 0, 2*sqrt(2)*x+2*z-sqrt(2) <= 0, 2*sqrt(2)*y-2*z-sqrt(2) <= 0, 2*sqrt(2)*y+2*z+sqrt(2) >= 0]), fraction)):
Plot(ps, color=grey, scaling=constrained, axes=normal, view=[(-1.4..1.4)$3]);  # The plot at the same scale along the axes

Here we show how to symbolically find the first 6 derivatives at  (0,1). If n> 6, then there are problems with the memory, because the expressions become too cumbersome.

p:= ln((1+x)*y) + exp(x^2*y^2) - x - cos(x);
fsolve(eval(p,y=1), x);
implicitdiff(p, y, x);
eval(implicitdiff(p, y, x), [x=0,y=1]);
seq(eval(implicitdiff(p, y, x$n), [x=0,y=1]), n=1..6);

Use  solve  command, which symbolically solves systems of polynomial equations in the form RootOf of a high degree polynomial, after that  allvalues  command and then apply  evalf  command.

Edit.

In your integral,  y  variable is the integration variable (dummy variable). A definite integral does not depend on the name of the variable of integration, so simply replace y under the integral sign with any symbol that is not in your worksheet. Here is a simple example in which the same error is encountered:

f := x->int(x/(x+sin(x)), x = 0 .. x^2); 
f(1/x);   # An error
F := x->int(s/(s+sin(s)), s = 0 .. x^2);  # Workaround
F(1/x);  # OK

Edit.

Example:

randpoly([x, y], dense, degree = 4);
select(t->degree(t,x)+degree(t,y)<2, %);  

Or even shorter:

select(t->degree(t, [x,y])<2, %);

Edit.