For the beautiful plot use in addition to the option  discont=true  also  the option  symbol=solidcircle  and some other options:

plot([floor(x), seq([i, t, t = min(0, i) .. max(0, i)], i = -3 .. 3)], x = -3 .. 3.9, thickness = [3, 1$7], color = red, linestyle = [1, 2$7], discont = true, symbol = solidcircle, symbolsize = 15);

N := 10; h := 1/N;

for i from 0 to N do x[i] := i*h; p[i] := evalf(cos(x[i])) end do:

f1 := [seq([x[i], p[i]], i = 0 .. N)];

plot([cos(x), f1], x = 0 .. 1, y = 0 .. 1, style=[line,point], color=[red,blue],symbolsize=15);

Instead of  from i = 1 to N do   you should write  

for i from 1 to N do  

or  

for i to N do  

(initial value  i  equals 1 by default)

For the rigorous proof of the absence of discontinuities in your example, you can use the identity

If the functions  f(x)  and  g(x)  are continuous,  this identity implies that  max(f(x), g(x))  is also continuous, because the function  abs  is continuous.

Obviously, the system  sys:= [eqn1-eqn2 = 0, eqn2-eqn3 = 0]  has the trivial solution  x=y=z=t  for any real  t . If I understand your goal, then you are trying to prove that this solution is the unique solution of the system. But this is not true, it is easy to see from the plot:

z := 1: 

sys:= [eqn1-eqn2 = 0, eqn2-eqn3 = 0]:

plots[implicitplot](sys, x = 0 .. 4, y = 0 .. 4, color = [red, blue], thickness = 2, gridrefine = 3);

My hypothesis  -  in typical case (excluding some individual values  t ) for each value of  z=t  the system  sys  has 4 solutions:[x=t, y=f(t), z=t],  [x=f(t), y=t, z=t],  [x=t, y=t, z=t], [x=g(t), y=g(t), z=g(t)]  where  f  and  g  are some functions.

This is probably a bug. In the procedures  Student[Calculus1]:-ShowSolution  and  Student[Calculus1]:-IntTutor  symbol  s  is automatically set to the empty list. Therefore, instead of the symbol  s , use any other symbol, for example,  a, b, r  and so on.

lista:=[[1,10],[2,20],[3,30]]:

x, y:=seq(map(t->t[i], lista), i=1..2);

                                            x, y := [1, 2, 3], [10, 20, 30]

restart;

z:=x->diff(y(x),x,x)+sin(y(x)):

u:=diff(z(x),x):

v:=diff(z(x),x,x):

w:=diff(y(x),x$3)+diff(y(x),x$2):

sys:=subs({y(x)=a, diff(y(x),x)=b, diff(y(x),x$2)=c, diff(y(x),x$3)=d, diff(y(x),x$4)=0.1},{z(x)=0, u=0, v=0, w=-0.1});

fsolve(sys);

assign(%);

sol:=dsolve({z(x)=0,y(0)=a,D(y)(0)=b}, numeric);

plots[odeplot](sol,[x,y(x)],x=0..3, thickness=2,view=[-0.5..3,-5..0.5]);

Addition: Only one solution is found. In fact, the original equation has four real solutions because the system  sys  has four real solutions  A:

restart;

sys := {c+sin(a) = 0, d+c = -0.1, d+cos(a)*b = 0, 0.1-sin(a)*b^2+cos(a)*c = 0};

A:=RealDomain[solve](sys);

Second addition: In fact, the set of solutions of  sys  will be infinite. Solutions are found only when  -Pi <a <Pi. If  {a, b, c, d}  is a solition then   {a+2*Pi*k, b, c, d} (k is an integer) is also a solution.

Plotting in classic Maple 12. To improve the quality, I added 2 options:

f:=x->max(x^2, sqrt(abs(x)));

plot(f, -3..3, thickness=2, numpoints=5000);

Addition: if you write

f:=proc(x)

if x^2>sqrt(abs(x)) then x^2 else sqrt(abs(x)) fi;

end proc;

then result is the same.

for arbitrary Matrix:

MatrixToString:=proc(A::Matrix)

local m, n;

m, n:=op(A)[1..2];

convert(cat(seq(op(["(", seq(op([convert(A[i,j], string), ","]), j=1..n), ")", "=="]), i=1..m-1), op(["(", seq(op([convert(A[m,j], string), ","]), j=1..n), ")"])), string);

end proc:

Example:

A:=Matrix(4,5, [$1..20]);

MatrixToString(A);

plots[implicitplot]  with  gridrefine  option also gives a nice plot:

plots[implicitplot]([x^3-4*x=y, y^3-3*y=x], x=-4..4, y=-4..4, color=[red,blue], thickness=2, gridrefine=3);

You do not need a graphics command. In the range  t=-0.05..0.05  the equation for each value  t  has a unique solution for  x , so  fsolve  command can be used:

data:=[seq([t, fsolve(t+sin(8*x^3*t)-x^3, x, maxsols=1)], t=-0.05..0.05, 0.001)];

plot(data, thickness=2);

P:=proc(n,x)

option remember;

simplify(((2*n-1)*x*P(n-1,x)-(n-1)*P(n-2,x))/n);

end proc:

P(0,x):=1:  P(1,x):=x:

P(6,x);

simplify(LegendreP(6,x));

testeq  command is not suitable for checking all the roots, you can use  eval  command.

Solution and check:

eqns := {x^3-4*x = y, y^3-4*y = x}:

vars := {x, y}:

solns := solve(eqns, vars, explicit):

solns1 := map(simplify, [solns]);  # Simplified solutions

seq(simplify(eval([lhs(eqns[1])-rhs(eqns[1]), lhs(eqns[2])-rhs(eqns[2])], solns1[i])), i = 1 .. 9); #Checking of all the solutions

You got the solutions in implicit form. Here is the complete list of the all solutions:

eqns := {x^3-4*x = y, y^3-4*y = x};
vars := {x, y};
solns := solve(eqns, vars);

map(allvalues, [solns]);