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

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