Your system can easily be solved exactly. In general case, it will have  2^3 = 8  solutions. Since the orthogonal matrix  RT has determinant  1, the matrix  mat  should also be orthogonal with determinant  1. For solution we need to select 3 independent elements of the orthogonal matrix  mat (I took  mat[1,2], mat[1,3], mat[2,3]). Here are all the solutions for the specific orthogonal matrix  mat  with full check of one solution:

restart;

Rx := Matrix(3,3, [1,0,0,0,cos(a),-sin(a),0,sin(a),cos(a)]):   

Ry := Matrix(3,3, [cos(b), 0, -sin(b), 0, 1, 0, sin(b), 0, cos(b)]):  # Corrected matrix

Rz := Matrix(3,3, [cos(c), -sin(c), 0, sin(c), cos(c), 0, 0, 0, 1]):

RT := Rx.Rz.Ry:

mat:=<1/3,2/3,-2/3; -2/3,2/3,1/3; 2/3,1/3,2/3>;  # The specific orthogonal matrix with det=1

assign(seq(seq(g[i,j]=(RT[i,j]=mat[i,j]), j=1..3), i=1..3));

Sol:=allvalues([solve({g[1,2], g[1,3], g[2,3]})]);  # Found solutions

eval([seq(seq(g[i,j], j=1..3), i=1..3)], Sol[1,1]);  # The check of  Sol[1,1]  for all equations of the system

Additions.  1) You can verify that  Sol[4,1]  is also a solution. The rest of the solutions are not solutions of the whole system.

2) The code was edited.

Your procedure works for me, but the same can be written a much shorter.

Your way:

dList:=proc(a::list, b::list)

    local res::list, i::integer, n::integer;

    description "add list elements and diff";

    n:=min(nops(a), nops(b)) ;

res:=[seq(0, i=1..n)];

for i from 1 by 1 to n do

     res[i]:=diff(a[i]+b[i],x)

end do;

res;

end proc:

list1:= [x^3 , 2*x , 3]:

list2:= [x^5, 2*x^4 , 2*x]:

dList(list1,list2);

                    [3*x^2+5*x^4, 2+8*x^3, 2]

The shorter method (without user's procedure):

list1:= [x^3 , 2*x , 3]:

list2:= [x^5, 2*x^4 , 2*x]:

map(diff, list1+list2, x);

                     [3*x^2+5*x^4, 2+8*x^3, 2]

 Maybe your error is caused by previous assignments. Make restart before your procedure.

Everything to be beautiful, I changed all the colors and ranges for the axes. Instead of  Student package I just wrote the simple equations for the tangent and normal:

f := x-> x^2-2*x:

x0 := 1.5:

T := (D(f))(x0)*(x-x0)+f(x0):  # The tangant

N := -(x-x0)/(D(f))(x0)+f(x0):  # The normal

P := plottools[point]([x0, f(x0)], symbol = circle, symbolsize = 15):  # The point

Curves := plot([f(x), T, N], x = -4 .. 5, -4 .. 5, color = [red, blue, green], thickness = 2, legend = [curve, tangent, normal]):

plots[display](Curves, P, scaling = constrained);

You will have more opportunities to vary the shape of the knot, if you use it's parametric equations and  plots[tubeplot]  command. You are setting the parameters  rho>1  and coprime integers  p  and  q .

Example:

restart;

r := cos(q*phi)+rho:  x := r*cos(p*phi):  y := r*sin(p*phi):  z := -sin(q*phi):  p := 3:  q := 7:  rho := 3:

plots[tubeplot]([x, y, z], phi = 0 .. 2*Pi, radius = 0.3, style = surface, numpoints = 3000, scaling = constrained, axes = none);

See Carl's idea using  viewpoint  option  here

In this thread are also discussed some other ways to solve the problem.

There is substantial error in procedure  Simpson  - instead  n  should be  n/2 . I made a few less significant corrections. Now everything is working correctly:

Simpson := proc(f, a, b, n::even)

local h, S1, S2, S, i;

h := (b-a)/n;

S1 := 0.;

for i from 0 to n/2-1 do

S1 := S1 + f(a + (2*i+1)*h);

end do;

S2 := 0.;

for i from 1 to n/2-1 do

S2 := S2 + f(a + (2*i)*h);

end do;

evalf(h/3 * ( f(a)+f(b) + 4*S1 + 2*S2 )); 

end proc:

Example of work of the procedure  Simpson  and independent calculation:

Digits := 5;                          

f := x -> 1/sqrt(39.24*x-44.65*(x*arccos(x)-sqrt(1-x^2)-13.88*(1-x^2)^1.5));

Simpson(f, 0, 1, 100);

evalf(Int(f(x), x=0..1));

is(diff(abs(x), x) = x/abs(x))  assuming x < 0;

is(diff(abs(x), x) = x/abs(x))  assuming x > 0;

                                  true

                                  true

Slightly more programmatically:

op(-1, eval(splcurve) assuming v < 4.2);

If you want to get a normal manual solution, you can use  Student[Calculus1][ShowSolution]  command.

Example:

Student[Calculus1][ShowSolution](Int(sin(x)*(x + sin(x)),x));

roll:=rand(1..6):

L:=[seq([seq(roll(), j=1..10)], i=1..30)];

NumbersOfThrees:=map(t->ListTools[Occurrences](3, t), L);

The number  pi  should be coded as  Pi .

All curves are in the same plot:

restart;

Eq:=ln(2*s*t+2*s*x*sqrt(t)/sqrt(Pi)+x^2)-2*s*(arctan((s+x*sqrt(Pi/t))/sqrt(2*Pi*s-s^2))-(1/2)*Pi)/sqrt(2*Pi*s-s^2) = 0;

S:=[0.01,0.005,0.003,0.002];

plots[implicitplot]([seq(eval(Eq,s=S[i]), i=1..nops(S))], t=0..5, x=-2..2, color=[red,blue,green,yellow], gridrefine=5);

It is easy to verify identity

is(sin(2*x)^2 - cos(8*x)^2 = - cos(10*x)*cos(6*x));

                                         true

Therefore, the equation splits into two equations, which are easy to solve

solve({x >= 0, `or`(cos(10*x) = 0, cos(6*x)+1/2 = 0), x <= (1/2)*Pi}, AllSolutions, explicit);

Consider the function which is defined by the equation and its symmetries:

restart;

f:=x->sin(2*x)^2-cos(8*x)^2-(1/2)*cos(10*x):

is(f(Pi-x)=f(Pi+x)) and is(f(Pi/2-x)=f(Pi/2+x));

                                        true

Therefore it is sufficient to solve this equation in the range  0 .. Pi/2

RootFinding[Analytic](f(x), x, re = 0 .. (1/2)*Pi, im = -1 .. 1);

identify([%]);

You can enter a variable (the counter), first assign it the value  0 , and then at each step of the loop is increased by 1 .

Example: find all Pythagorean triangles with hypotenuse not exceeding 100 . Solution in the for loop:

restart;

k:=0: # The counter for the steps of the loop

n:=0: # The counter for solutions

for c from 2 to 100 do

for a from 1 to c-1 do

k:=k+1:

b:=sqrt(c^2-a^2):

if type(b,integer) and b>a then n:=n+1; L[n]:=[a,b,c] fi;

od:

od:

L:=convert(L, list):

k;  n;  L;

Addition:  The example of a purely illustrative nature. There are far more effective methods of generating Pythagorean triples.

restart;

 L:=[seq(seq(x^i*y^j, j=0..3-i), i=0..3)];

 add(a[k]*L[k+1], k=0..nops(L)-1);