I am struggling with a definition of a functional surface z=z(x,y) where z is positive in the interior of a triangle, zero on its boundary, and zero outside of it everywhere on a rectangle for (x,y) which contains the triangle. All necessary formulae can be provided in if/elif/else/return, etc. statements. Yet, I get errors when using def and return (I have never worked with simpy and would prefer to avoid using lambdas and stay as close to the mathematical expressions as possible). I think that if someone helps with explaining how to handle a simple definition of a 1-variate function of x=var('x'), I'll figure what I am not doing right also in the 2-variate case x, y = var('x,y'). So here is a very simple example:

f has definition domain -2<=x<=3 and is defined with different formulae on different parts of its domain, e.g.,

f(x)=x+1 when x is between -2 and -1,

f(x)=0 when x is between -1 and 0,

f(x)=x^2 when x is between 0 and 1,

f(x)=x^3 when x is between 1 and 3,

(naturally, without overlapping at the boundaries, although this function is continuous there, so no problems arise there -- I skipped the use of inequalities because the preview program got confused).

After the definition I would like to use f(x) while specifying for x the WHOLE range [-2,+3]. (For example, plot the graph of y=f(x) for all x between -2 and +3.) Can someone please provide a simple working code, say, using def and conditional statements but not resorting to tuples, classes (ok, if there is absolutely not a simple way, lambdas can be used as a last resort :) )? Or maybe this question has been answered earlier? Please help :)

Following the fairly exhaustive answers of rburing and tmonteil below to the 1-variate version of my question, I now proceed to upgrading it for 2-variate functions with a maximally simplified concrete problem which, if solved, can immediately be upgraded also to the 3-variate case. (In the formulation, I shall use latex notation.)

Consider the square domain (x,y) \in A=[-2,2]x[-2,2] and a subset of its interior: the 'canonical' triangle B with vertices P_0=O=(0,0), P_1=(1,0) and P_2=(0,1). The reason for this 'canonical' choice is that the cartesian coordinates x and y coincide with two of the three barycentric coordinates of an arbitrary point P in the plane Oxy:
P=(x,y)=(1-x-y).P_0+x.P_1+y.P_2; moreover, if P is in B or on its boundary, then the combination is CONVEX; (moreover, if P is on the interior of an edge of the boundary of B, then 1 of the three coefficients is 0, and if it is coincides with one of the vertices, then 2 of these coefficients are 0). Now I proceed to define mathematically the needed functional surface z=z(x,y), as follows:

Consider first a point P_3=(x_3, y_3) which lies in the interior of B, i.e. all the three coefficients 1-x-y, x and y are strictly positive and strictly less than 1.

if 1-x-y leq 0 or x leq 0 or y leq 0 (this is exactly the case when P=(x,y) is outside of B or lies on its boundary), then z(x,y)=0;

elsewhere on A (i.e., exactly when P lies in the interior of B), define
z(x,y)=exp{-[(x-x_3)^2+(y-y_3)^2]/[(1-x-y)^2(x^2)(y^2)]}.

Now z(x,y) is defined on the whole A and has several remarkable properties there. As a non-exhaustive example:

it is analytic everywhere on A except the boundary of B where it is still C^\infty - smooth but not analytic (the Taylor expansion series is not convergent there, but is convergent everywhere else on A).

it takes strictly positive values less than 1 (good as coefficient in a convex combination!) and reaches its global maximum equal to 1 exctly and only in the prescribed point P_3.

To define this type of functional surface on an arbitrary triangle in Oxy or as a parametric surface in 3D one now uses a standard approach (used in finite/boundary element theory, geometric design, etc, etc.) via applying a non-singular affine transformation. The technique is easy to extend to tetrahedra in 3D.

To explain why this function has numerous prospective valuable applications is a very long story. To cut it short, I promise that if you help solve this particular problem, you will soon read about some of its new applications in the newspapers... :) :) :)

Summarizing:

the primary, maximally simplified, task I ask to be solved here is to provide SageMath (Python 3) code to use parametric_plot3d to plot the parametric surface X(u,v)=u, Y(u,v)=v, Z(u,v)=z(u,v) for (u,v) in the parametric domain A given above and to view it in Ubuntu using threejs and/or jmol. (Forget about all the extra options, Windows limitations, etc).

Can you help, please?