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.

Simple procedure Cramer returns the main determinant of the system  det[A] , all the auxiliary determinants  det[i]  and the values of all unknowns calculated according to Cramer's rule:

Cramer:=proc(Sys::{set(equation), list(equation)}, Var::list(name))
local n, A, V, det, S;
uses LinearAlgebra;
n:=nops(Sys);
A, V := GenerateMatrix(Sys, Var);
det[A]:=Determinant(A);
S:=seq(Determinant(<A[..,1..i-1] | V | A[..,i+1..n]>), i=1..n);
print(det__A=det[A], seq(det[i]=S[i], i=1..n));
seq(Var[i]=S[i]/det[A], i=1..n);
end proc:

Example of use:

A:=LinearAlgebra:-RandomMatrix(4, generator=-9..9):
V:=LinearAlgebra:-RandomVector(4, generator=-9..9):
Sys:=LinearAlgebra:-GenerateEquations(A, [seq(x[i], i=1..4)], V);
Cramer(Sys, [seq(x[i], i=1..4)]);
    
 

restart;
with(DEtools): 
DE1 := diff(y(x),x) = diff(y(x),x$3)+3*diff(y(x),x$2)+4*diff(y(x),x)+12*y(x);
DEplot(DE1, y(x), x = 0..5, [[y(0)=-3, D(y)(0)=0, (D@@2)(y)(0)=0]]);

Your equation only Maple 2017 will be able to solve:

pde := (diff(u(x, y), x, x))+diff(u(x, y), y, y) = 0;
 bc[1]:= D[2](u)(x,0) = 0;
 bc[2]:= u(x,1) = x^2-x;
 bc[3]:= u(0,y)=0;
 bc[4]:= u(1,y)=0;
 pdsolve({pde, bc[1], bc[2],bc[3],bc[4]}, HINT= X(x)*Y(y));

HINT  option can be omitted.

X:=<1,2,3,4>:
Y:=<5,6,7,8>:
x.X+y.Y;

In this example  X  and  Y  are vectors. For matrices everything is the same.

Edit.

Try  ContoursWithLabels  procedure  from  here .

Your example:

ContoursWithLabels(-3.392318438*exp(-4.264628895*x)*sin(6.285714286*y), -1/2 .. 1/2, -1/2 .. 3/2, {seq(-12 .. 12, 3)}, [color = black, axes = box, size=[1000,600]], Coloring = [colorstyle = HUE, colorscheme = ["Cyan", "Red"], style = surface]);

If you do not want coloring, then just remove  Coloring  option from the code.

With the built-in command  plots:-contourplot , do so

plots:-contourplot(-3.392318438*exp(-4.264628895*x)*sin(6.285714286*y), x = -1/2 .. 1/2, y = -1/2 .. 3/2, contours=[seq(-12..12,3)], size=[800,500], numpoints=10000);

Edit.

Execute:

op~(Sol);