It makes no sense to look for what does not exist, because no global minimum and global maximum. A point  (x,y,z)  lies on the sphere of radius  sqrt(5)  and the plane  x+y+z+1=0   intersects the sphere. Therefore, the denominator of the fraction  96/(x+y+z+1)  can be arbitrarily close to 0 and may be positive or negative.

This can be seen clearly in the plot:

phi:=Pi/4:
A:=subs(x=sqrt(5)*sin(theta)*cos(phi), y=sqrt(5)*sin(theta)*sin(phi), z=sqrt(5)*cos(theta), 1/2*x^2*y^2 + y^2*z^2 + z^2 *x^2 + 96/(x + y + z + 1)):
plot(A, theta=0..Pi, view=[0..Pi, -1000..1000], thickness=2, discont=true);

Should be

u:=(x,t)->exp(-Pi^2*t)*cos(Pi*x):

Your problem can be solved by the usual brute force. The program returns a list of all suitable  [a, b​​, c]  and the corresponding sequences  t[k], k=-1..20

t[-1]:=0: t[0]:=1: L:=[]: T:=[]:

for a to 20 do

for b to 20 do

for c to 20 do

for n from 0 to 19 do

t[n+1]:=((a*n^2+a*n+b)*t[n]+c*(n^2)*t[n-1])/(n+1)^2:

od:

if convert([seq(type(t[k], integer), k=1..20)], `and`) then L:=[op(L), [a,b,c]]:

T:=[op(T), [seq(t[k], k=-1..20)]] fi:

od: od: od:

L; T;

Two errors in your code. Instead of  (H[1]-x2)*(x2-x1) +(H[2]-y2)*(y2-y1)=0 and (x2-x1)*(y3-y1) +(x3-x1)*(y2-y1) <>0  should be  (H[1]-x3)*(x2-x1) +(H[2]-y3)*(y2-y1)=0 and (x2-x1)*(y3-y1) -(x3-x1)*(y2-y1) <>0 

Let the vertices of the triangle will be  [x1,y1], [x2,y2], [x3,y3] . Firstly we find the coordinates of orthocenter:

restart:

sol:=solve({(x3-x2)*(x-x1)+(y3-y2)*(y-y1)=0, (x2-x1)*(x-x3)+(y2-y1)*(y-y3)=0}, {x,y});

assign(sol):

Next, we find the solution of the original problem. Note that can confine ourselves 4-loops instead of 6-loops, if we use the formula for the centroid. So we can increase the range of  search. Also, instead of lists we use sets, so exclude the same triangles (which are obtained by permutations of the vertices):

L:={}:

k:=50:

for x1 from -k to k do

for x2 from -k to k do

for y1 from -k to k do

for y2 from -k to k do

x3:=3-x1-x2: y3:=3-y1-y2:

if x2*y3-x1*y3+x1*y2-y2*x3+y1*x3-y1*x2<>0 and x=3 and y=3 then

L:={op(L), {[x1,y1], [x2,y2], [x3,y3]}}: fi:

od: od: od: od:

nops(L); L;

If the centroid of a triangle coincides with the orthocenter, then the triangle is equilateral. It is easy to prove that in 2D coordinates of all the vertices of this triangle can not be integers.

int(rho*(rho+1)*sin(phi), [rho = 0 .. R, phi = 0 .. Pi, theta = 0 .. 2*Pi]);

Tridiag:=proc(Diag1, Diag, Diag2, n)

local m, A, i, j;

m:=nops(Diag);

A:=Matrix(m, n);

for i to m do

for j to n do

if i=j then A[i,j]:=Diag[i] else

if i=j+1 then A[i,j]:=Diag1[j] else

if j=i+1 then A[i,j]:=Diag2[i] else

A[i,j]:=0; fi; fi; fi;

od; od;

A;

end proc:

Example:

Tridiag([1,2,3,4,5], [2,3,4,5,6,7], [3,4,5,6,7,8], 8);

We assume that the cylinder of a revolver has 6 chambers. The procedure returns a structure of each trial and vector of   times the trigger was pulled.

RR:=proc(N)  # N is the number of trials

local roll, i, a, L, j, b, M, Q, A;

roll:=rand(1..6);

for i to N do

a:=roll(); L:=[[a]];  # a is the chamber number in which the bullet

for j do

b:=roll(); if b=a then L:=[op(L), b]; break else  L:=[op(L), b]; fi;

od;

M[i]:=L;

od;

Q:=[seq(M[i], i=1..N)];

A:=Vector([seq(nops(Q[i])-1, i=1..nops(Q))]);

print(op(Q));

print(A);

end proc:

Example:

RR(10);

SolvingSystem:=proc(system, n)
 if n=1 then method s1 fi;
 if n=2 then method s2 fi;
 if n=3 then metod s3 fi;
 result;
end proc:

add(add(int00[i,j], j=1..3), i=1..3);

                           0

You have to write  exp(-10*x)  instead of  exp^(-10*x) 

Should be:

l4:=animate(textplot3d, [[0, 0, 0.9*tend,'direction'], color=red, align={BELOW}, font=[TIMES,ROMAN,14]], tend=0..50, frames=150):

 Coordinates of the end of a vector are the coordinates of the basis  plus the coordinates of the vector.

See  help on commands  ?plots[matrixplot]  and  ?table

code1.mw

Your error: in the last lines should be the sign  :=  instead of  =

Let  x  is serial number 5, y  is serial number 6, z  is serial number 7.

L:=[]:
for x from 1 to 2 do
for y from 0 to 9 do
for z from 0 to 9 do
 if z=5-y then L:=[op(L), [x,y,z]]: fi:
od: od: od:
L; nops(L); combinat[numbperm]([9,8,7,7])*%;
 

 [[1, 0, 5], [1, 1, 4], [1, 2, 3], [1, 3, 2], [1, 4, 1], [1, 5, 0],

        [2, 0, 5], [2, 1, 4], [2, 2, 3], [2, 3, 2], [2, 4, 1],

        [2, 5, 0]]

                                  12

                                 144