Obviously, the problem reduces to the calculation  cos(z)  and  sin(z) , if  z=1/2*arccos(1/3) . Using Euler's formula  exp(I*z)=cos(z)+I*sin(z), we obtain

expr:=expand(sin(-(1/6)*Pi+(1/2)*arccos(1/3)));

z:=1/2*arccos(1/3):

evalc(exp(I*z));  # Euler's formula in Maple

eval(expr,[cos(z)=Re(%),sin(z)=Im(%)]);

The easiest way to draw graphs of numerical solutions of differential equations - is to use  plots:-odeplot  command.

Example with 2 functions:

Sol:=dsolve({diff(y[1](x),x)=sin(y[2](x))+2*cos(y[1](x)), diff(y[2](x),x)=2*cos(y[2](x))+3*sin(y[1](x)), y[1](0)=0, y[2](0)=1}, numeric):

plots:-odeplot(Sol, [[x,y[1](x)], [x,y[2](x)]], x=0..4, color=[red,blue], thickness=2, view=[0..4, 0..4]);

Addition. If you still want to use  plot  command, you can write as in example: 

y[1]:=x->rhs(Sol(x)[2]):

plot(y[1], 0..4);

Your contour - it is an ellipse. To maintain the proportions along the axes Ox and Oy is necessary to take the same ranges. For clarity, the cutting plane  z=160   is also shown:

A := plots:-contourplot3d(4*x^2+9*y^2, x = -8.5 .. 8.5, y = -8.5 .. 8.5, contours = [160], color = red, thickness = 4, linestyle = solid):

B := plot3d(4*x^2+9*y^2, x = -8.5 .. 8.5, y = -8.5 .. 8.5, style = surface, view = -15 .. 190):

C:=plot3d(160, x=-8..8, y=-6..6, transparency=0.5):

plots:-display(A, B, C, axes = normal, orientation = [55, 55], lightmodel = light4);

Below we prove this using Maple for calculations.
First, we prove that the number  n  can not have more than 21 digits.

The equation  9^20*(n+1)=10^n  has not more than 2 roots as the intersection of a straight line with a convex curve. Find these roots and the plot (for clarity) :

Digits:=25:

fsolve(9^20*(n+1) = 10^n, n = -2 .. 0);

fsolve(9^20*(n+1) = 10^n, n = 0 .. infinity);

plot([9^20*(n+1), 10^n], n = -2 .. 22, 0..4*10^20, thickness=2);

Next, check that the number of 199 ... 9 (20 nines) satisfies the required property:

is(1+9^20*20>=10^20+10^20-1);

                           true

It remains to verify that any 21-digit number greater than 199 ... 9, does not satisfy this property. There is no need to check all the numbers up to 10^21-1. The basic idea - the sum  20th powers of the digits from the number  199 ... 9  to the number  99 ... 9  rises very slowly (just 1.05 times), so to check the property sufficiently in the small range. Changing 4 digits and doing the simple estimates we get the desired result:

restart;

evalf(9^20*21/(1+9^20*20));

evalf((1+9^20*20)/(10^20+10^20-1));

%*%%;

(2*10^20-1)*%;  # The maximum number up to which we should to check

s:=1: 

for a from 0 to 5 do

for b from 0 to 5 do

for c from 0 to 3 do

for d from 0 to 2 do

N1:=2*10^20+a*10^19+b*10^18+c*10^17+d*10^16:

N2:=2*10^20+a*10^19+b*10^18+c*10^17+d*10^16+10^16-1:

if  2^20+a^20+b^20+c^20+d^20+16*9^20 < N1 then s:=s*1 else s:=s*0 fi;

od: od: od: od:  

s;

with(numapprox):

pade(cos(x), x, [4,4]);

plot([cos(x), %], x=0..2*Pi, -1.2..1.5, color=[red,blue]);

We see that in the range  x=0..3  approximation is very good - the curves merge.

Obviously, this number cannot have 5 or more digits because 9^3*5<10000. 

Simple enumeration solves the problem:

restart;

for k to 9999 do

L:=convert(k,base,10);

s:=add(L[i]^3, i=1..nops(L));

if s>=k then n:=k fi;

od:

n;

                         1999

simplify(abs(x)-sqrt(x^2))  assuming x::real;

                                    0

This issue differs from the previous in that the right-hand side of the equation is much smaller. Therefore, the possibility appeared to estimate the  number of solutions. It turned out that solutions more than 100000000. Furthermore, among the solutions, we find the solution with minimum norm.

Finding the number of solutions:

restart;

eq:=30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

isolve(eq);

ineq:=5*_Z1+10*_Z2+25*_Z3+35*_Z4+25*_Z5+40*_Z6+5*_Z7+10*_Z8<640;

L:=coeffs(lhs(ineq));

k:=0:

for _Z1 to (640-`+`(L[2..8]))/5 do

for _Z2 to (640-5*_Z1-`+`(L[3..8]))/10 do

for _Z3 to (640-5*_Z1-10*_Z2-`+`(L[4..8]))/25 do

for _Z4 to (640-5*_Z1-10*_Z2-25*_Z3-`+`(L[5..8]))/35 do

for _Z5 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-`+`(L[6..8]))/25 do

for _Z6 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-25*_Z5-`+`(L[7..8]))/40 do

for _Z7 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-25*_Z5-40*_Z6-L[8])/5 do

for _Z8 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-25*_Z5-40*_Z6-5*_Z7)/10 do

if ineq then k:=k+1 fi;

od: od: od: od: od: od: od: od:

k;  # The  number of solutions for positive integer parameters  _Z1, _Z2, ..., _Z8

Next we find a solution with minimum norm. To do this, first we find the optimal solution in two ways without regard integer, and then varying the closest integer solutions solves the problem:

restart;

eq:=30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

F:=a^2+b^2+c^2+d^2+e^2+f^2+g^2+h^2+i^2:

Optimization[QPSolve](F,{eq}, assume=nonnegative);

extrema(F, {eq}, indets(F),'s'): evalf(convert(op(s),list));

map(t->[round(rhs(t)-4),round(rhs(t)+4)], %);

Dist:=infinity:

for a from 1 to 6 do

for b from 1 to 9 do

for c from 3 to 11 do

for d from 1 to 9 do

for e from 12 to 20 do

for f from 6 to 14 do

for g from 1 to 5 do

for h from 1 to 5 do

for i from 3 to 11 do

if eq then dist:=F;

if dist<Dist then Dist:=dist; Sol:=[a,b,c,d,e,f,g,h,i]  fi; fi;

od: od: od: od: od: od: od: od: od:

['a'=Sol[1],'b'=Sol[2],'c'=Sol[3],'d'=Sol[4],'e'=Sol[5],'f'=Sol[6],'g'=Sol[7],'h'=Sol[8],'i'=Sol[9]];  # The solution with minimal norm

Edit. Comments on the first  code was corrected. In fact, this code finds not all solutions but only part of the solutions for positive values of parameters. The total number solutions is  more than 106141495 solutions.

Example:

plot([x,x^2,x^3], x=-1..1, thickness=2, legend=[x,x^2,x^3]);

for your way:

unapply(diff(f(x), x), x);

                  x->2*x

Of course, the use of the operator  D  is more preferably.

M:=8:  K:=6:  K1:=K-2*K2:

Sum(M!/(M-K1-K2)! * K!/(K1! * K2)! * 1/M^K * 1/2^K2, K2=0..K/2);

value(%);

I do not understand why you need something to simplify. You can easily to define  this matrix as done below and calculate its determinant with any precision for any  k  and  alpha :

Example:

L:=[J,Y,`I`,K]:

A:=Matrix(4,(i,j)->`if`(i=1 or i=2,cat(Bessel,L[j])(i+1,2*k*sqrt(alpha)/(alpha-1)),cat(Bessel,L[j])(i+1,2*k/(alpha-1)))):

F:=unapply(Matrix(4,{seq(seq(`if`(j=3 and (i=2 or i=4),(i,j)=-A[i,j], (i,j)=A[i,j]),j=1..4),i=1..4)}), k,alpha):

F(k,alpha);

evalf[20](LinearAlgebra[Determinant](F(2,3)));

 Edit. The code and the example was corrected (at first I missed minus signs)

I corrected your syntax as I understood it:

f := x->piecewise(abs(x)<1,1-x,0);

U := (x,t)->2/3*f(x+t)+1/3*f(x-2*t)+1/(3*Pi)*sin(Pi*(x+t))+1/(3*Pi)*sin(Pi*(x-2*t));

plot3d(U(x,t), x=0..3, t=0..3, axes=normal, numpoints=10000);

The original problem has an obvious meaning, if one or both parameters  H1  and  H2  are 0 or negative, but the proposed solution does not work for these cases. Here is the solution in the form of a procedure that works for all cases. The procedure is very simple, it is based on direct calculation of the volume by a single integration over the cross sectional area (each a cross-section - this is the usual circle segment).

restart;

V:=proc(R,H1,H2)

local r, H;

if is(H1^2+H2^2>R^2) then error "Should be H1^2+H2^2<=R^2" fi;

r:=int((R^2-x^2)*arccos(H1/sqrt(R^2-x^2))-H1*sqrt(R^2-x^2-H1^2), x=H2..sqrt(R^2-H1^2));

if is(H1>=0) then return simplify(r) else

H:=R-sqrt(R^2-H1^2);

simplify(r+Pi*H^2*(R-H/3)) fi;

end proc:

Examples of use:

V(1,1/2,1/2);

 evalf(%);

V(1,0,1/2);

V(1,1/2,0); 

V(1,-1/2,1/2);

 evalf(%);

V(1,1/2,-1/2);

 evalf(%); 

V(1,0,0); 

V(1,-1/2,-1/2);

 evalf(%);

V(1,-1/sqrt(2),-1/sqrt(2));

 evalf(%);

Here are solutions the same examples by Rouben's formula:

Vwedge :=(R,a,b)->-( (1/3)*Pi*R^3-(1/6)*Pi*(2*R+a)*(R-a)^2-(1/6)*Pi*(2*R+b)*(R-b)^2-(1/3)*a*(3*R^2-a^2)*arcsin(b/sqrt(R^2-a^2))-(1/3)*b*(3*R^2-b^2)*arcsin(a/sqrt(R^2-b^2))-(2/3)*a*b*sqrt(R^2-a^2-b^2)+(1/3)*R^3*arctan((R^2+R*b-a^2)/(sqrt(R^2-a^2-b^2)*a))-(1/3)*R^3*arctan((R^2-R*b-a^2)/(sqrt(R^2-a^2-b^2)*a))):

Vwedge(1,1/2,1/2);

 evalf(%);

Vwedge(1,0,1/2);

Vwedge(1,1/2,0);

Vwedge(1,-1/2,1/2);

 evalf(%);

Vwedge(1,1/2,-1/2);

 evalf(%);

Vwedge(1,0,0);

Vwedge(1,-1/2,-1/2);

 evalf(%);

Vwedge(1,-1/sqrt(2),-1/sqrt(2));

Edit.  The procedure  V  worked incorrectly if  H1<0 . So the procedure and all the examples was corrected. Now all right. The number of examples was increased. 

plot3d([sin(u)*(7+cos(u/3-2*v)+2*cos(u/3+v)), cos(u)*(7+cos(u/3-2*v)+2*cos(u/3+v)),

sin(u/3-2*v)+2*sin(u/3+v)], u=-Pi..Pi, v=-Pi..Pi,  axes=normal, orientation=[30,30], lightmodel=light4);

Edited: the extra commas was removed.