See corrected file   command_syntax_new.mw

Edit.

If your code is left as it is, then for large values  j  huge symbolic expressions will accumulate, which very slows down the work. So I instead of  rand()  wrote  rand(0...evalf(2*Pi)) , which gives a random angle in the range  0..2*Pi. I also took  c=1. All this greatly accelerated the work of the code. For visualization, I took the first 2 points and 100 steps. With  randomize()  command each time we run this code, we can see a different picture.

 restart;
M:=10: N:=1000: c:=1: r:=rand(0...evalf(2*Pi)):
randomize():
for i to M do
X[i, 0] := 0;
Y[i, 0] := 0;
X[i, 1] := 1;
Y[i, 1] := 0;
for j from 2 to N do
Vinkel := r();
X[i, j] := X[i, j-1]+cos(Vinkel);
Y[i, j]:= Y[i, j-1]+sin(Vinkel)
end do:
end do:
plot([[seq([X[1,j],Y[1,j]], j=0..100)], [seq([X[2,j],Y[2,j]], j=0..100)]], color=[red,blue], scaling=constrained, axes=box);  # Plotting for two points
j:='j':
plots:-animate(plot,[[['seq'([X[1,j],Y[1,j]], j=0..round(s))], ['seq'([X[2,j],Y[2,j]], j=0..round(s))]], color=[red,blue], scaling=constrained, axes=box], s=0..100, frames=200, paraminfo=false);  # Animation for two points 

The plot for  j=0..1000 :

 plot([[seq([X[1,j],Y[1,j]], j=0..1000)], [seq([X[2,j],Y[2,j]], j=0..1000)]], color=[red,blue], scaling=constrained, axes=box);

Edit.

Examples:

sol1 := solve({x+y+z=1, x-y+3*z=7});
sol2 := solve({x+s*y+z=1, x-y+3*z=7*t}, {x,y,z});
Params1:=indets(map(rhs, sol1), name);
Params2:=indets(map(rhs, sol2), name);

Here's the first option as I first understood the problem:

restart;
M0:=100: R3:=rand(1..3):
X[1] := 1: Y[1] := 0:
randomize():
for i from 2 to M0 do
r := R3();
if r = 1 then X[i] := X[i-1]; Y[i] := Y[i-1]+1;
elif r = 2 then X[i] := X[i-1]+1; Y[i] := Y[i-1];
else X[i] := X[i-1]-1; Y[i] := Y[i-1];
end if;
end do:

A:=plots:-animate(plot,[['seq'([X[j], Y[j]],  j = 1 .. round(n))], linestyle=2, color=blue, scaling = CONSTRAINED], n=1..M0, frames=180,paraminfo=false):
B:=plots:-animate(plots:-display,['plottools:-disk'([X[round(n)], Y[round(n)]],0.2, color=red)], n=1..M0, frames=180, paraminfo=false):
plots:-display(A, B, axes=box);

        

Here is another option (probably more correct, see tomleslie's answer and my comment to it):

restart;
M0:=100: R3:=rand(1..3):
X[1] := 0: Y[1] := 0: X[2] := 1: Y[2] := 0:
randomize():
for i from 3 to M0 do
r := R3();
if r=1 then X[i]:=2*X[i-1]-X[i-2]; Y[i]:=2*Y[i-1]-Y[i-2];
elif r=2 then X[i]:=X[i-1]+Y[i-1]-Y[i-2]; Y[i]:=Y[i-1]-X[i-1]+X[i-2];
else X[i]:=X[i-1]-Y[i-1]+Y[i-2]; Y[i]:=Y[i-1]+X[i-1]-X[i-2];
end if;
end do:

A:=plots:-animate(plot,[['seq'([X[j], Y[j]],  j = 1 .. round(n))], linestyle=2, color=blue], n=1..M0, frames=180,paraminfo=false):
B:=plots:-animate(plots:-display,['plottools:-disk'([X[round(n)], Y[round(n)]],0.2, color=red)], n=1..M0, frames=180, paraminfo=false):
plots:-display(A, B, axes=box, view=[min(seq(X[i],i=1..M0))-1..max(seq(X[i],i=1..M0))+1,min(seq(Y[i],i=1..M0))-1..max(seq(Y[i],i=1..M0))+1], scaling = constrained);

     

Edit.

Random_walk.mw

It's true for any real  n :

expr:=exp(-n*ln(2*Pi));
simplify(expr, exp) assuming n::real;

# Or
RealDomain:-simplify(expr, exp);

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