It is interesting to compare  two presented approaches to the solution of the problem. If you want to portray the whole cone, Markiyan's approach is preferable because it is more simple and does not require any knowledge of the equations of the cone. But if you want to construct only part of the cone, then I do not know how to do it using primitives. My method is more versatile because it allows you to do this easily. For example, if you need to cut a quarter of the cone , it is enough just to change 2 numbers in the code:

R:=1:  H:=2:

A:=plot3d([r*cos(phi), r*sin(phi), H-H/R*r], r=0..R, phi=Pi/2..2*Pi):  # The lateral surface

B:=plot3d([r*cos(phi), r*sin(phi), 0], r=0..R, phi=Pi/2..2*Pi):  # The bottom

plots[display](A, B, axes=normal, scaling=constrained);

@lcrpn   sin(Pi/16)=convert(sqrt((1-cos(2*Pi/16))/2), radical);

@lcrpn   The result is expressed in radicals (Classic Maple 12):

Rule:=a::realcons*sin(x)+b::realcons*cos(x)=sqrt(a^2+b^2)*sin(x+arctan(b,a)):

simplify(applyrule(Rule, collect(sin(x)+convert(expand(cos(x+Pi/5)), radical), sin(x))));

Addition:

In fact, this expression can be simplified considerably more (easier to do it manually and check in Maple):

A:=sin(x)+cos(x+Pi/5)=2*sin(3*Pi/20)*sin(x+7*Pi/20);

B:=lhs(A)=convert(2*sin(3*Pi/20),radical)*sin(x+7*Pi/20);  #  In radicals

is(A);  # Checking of equality in Maple

@Gauss  Your designation  NullSpace^-1(Y)  does not make sense, since it is neither a matrix or operator, but  a subspace, so just write  NullSpace(Y) .

If  Y  is a known matrix, for example, such as  Y=<< 1 | -2 | 1>, <-6 | 3 | 0>, <4 | -5 | 2 >>, then to find NullSpace(Y) necessary to solve the homogeneous system of linear equations:

restart;

Y:=<<1 | -2 | 1> , <-6 | 3 | 0> , <4 | -5 | 2>>;

LinearAlgebra[LinearSolve](Y, <0,0,0>);

We see that the NullSpace(Y) is one-dimensional subspace of  R^3  spanned by the vector  <1, 2, 3>

@Markiyan Hirnyk  I mean, if  topic-starter wants to obtain a unique solution, but not an infinite set of solutions, expressed in terms of the extra unknowns.

@Markiyan Hirnyk   Yes, of course

convert(exp(x), Sum, x=1, dummy=n);

@Gauss   but for a given matrix  A  the solution is unique, and your  NewMatrix3  is contrary to found  solution.

@PatrickT   I deleted the extra bracket.

Picture of output loaded in another way.

@Gauss   On my machine, the code is executed in a few seconds (in Classic M 12).  Run  restart  command  first.

@toandhsp  No, it can not. AreIsometric  procedure works only with the sets in the Euclidean plane, rather than in space.

@Alejandro Jakubi  Thanks for the explanation of the reason of the error.

@taro yamada  My method is only suitable for Classic Worksheet. If you write

A:=x/2+y/2:

normal(A/a);

then you see    

And then you just replace  symbol  a  by  the empty symbol  `` and get the desired output:

In the previous embodiment, the diagonal entries are always  0 , and should be random integers from  -9  to  9 . New version is free from this drawback:

roll:=rand(-9 .. 9):
Matrix(4, {seq((i, i)=roll(), i=1..4), seq(seq((i, j)=roll()+I*roll(), j=i+1..4), i=1..3)}, shape=hermitian);

Factoring  procedure takes 1/2 sec for this calculation, and the improvements made by Joe Riel and Carl Love even less time. See  http://www.mapleprimes.com/posts/141668-Partitions-Of-A-Natural-Number-Into-Factors

ts:=time():

Factoring(9!):

nops(%);

%%[1..20];

time()-ts;

I have already answered a similar question about the matrix  A . About the matrix  B  all is still easier. There is an inaccuracy in determining  B. Apparently instead of  -1/(2*M*(M-2))   should be -1/(M*(M-2)) , otherwise a contradiction with special cases.