#From cheung@math.lsa.umich.edu Fri Sep 10 12:57:58 1993
#To: carlos


mvcal[approxint2d] := proc(f,xexpress, yexpress ,n ,m)
           local a,b,c,d, dx,i,j, dy;
        a := op(1, op(2, xexpress));
        b := op(2, op(2, xexpress));
        c := op(1, op(2, yexpress));
        d := op(2, op(2, yexpress));
                         dx := evalf((b-a)/n);  dy := evalf((d-c)/m);
          sum(sum(subs(x=a+i*dx,y=c+j*dy,f),i=1..n),j=1..m)*dx*dy;
end: 


mvcal[arrows] := proc(vector1, vector2) 
local v1, v2, u1, u2 ;
    v1 :=  vector1[1]; v2 := vector1[2]; u1 := vector2[1] ; u2:= vector2[2] ;
    plot({ [0,0, v1, v2], [v1, v2, (4/5*v1-v2/10), (4/5*v2+v1/10)], [v1, v2,
(4/5*v1+v2/10), (4/5*v2-v1/10)],[0, 0, u1, u2], [u1,
u2,(4/5*u1-u2/10),(4/5*u2+u1/10)], [u1, u2,(4/5*u1+u2/10), (4/5*u2-u1/10)] }) ;
end :

mvcal[arrows3d] := proc( vector1, vector2)
local u1 , u2, u3,   v1, v2, v3;
u1 := vector1[1] ; u2 := vector1[2]; u3 := vector1[3]; v1 := vector2[1];
v2 := vector2[2]; v3 := vector2[3];  
plot3d({[t*u1, t*u2, t*u3], [t*v1, t*v2, t*v3], [u1-s/2*u1, u2-s*u2, u3*(1-s)],
[u1-s*u1, u2-s*u2/2, u3*(1-s/2)] , [v1-s/2*v1, v2-s*v2, v3-s*v3], [v1-s*v1,
v2-s*v2/2, v3-s*v3/2]}, t=0..1, s=0..0.3, grid=[2,2]); end :


mvcal[blockapp] := proc( f, xexpress, yexpress)
local x1, x2, y1,y2, i ,j , a, b ,c, d:
a := op(1, op(2, xexpress));
 b := op(2, op(2, xexpress));
 c := op(1, op(2, yexpress));
  d := op(2, op(2, yexpress));
blocks := {} :
for i from 1 to 4 do
x1 := a + (i-1)*(b-a)/4:
x2 := a + i*(b-a)/4:
for j from 1 to 6 do
y1 := c + (j-1)*(d-c)/6:
y2 := c + j*(d-c)/6:
blocks := blocks union {[x1, s*y1+ (1-s)*y2, t*subs(x=x2,y=y2, f)], [x2,
s*y1+(1-s)*y2, t*subs(x=x2, y=y2, f)], [s*x1+ (1-s)*x2, y1,
t*subs(x=x2,y=y2, f)], [ s*x1+ (1-s)*x2, y2, t*subs(x=x2,y=y2, f)],
[s*x1+ (1-s)*x2, t*y1+(1-t)*y2, subs(x=x2,y=y2, f)]} :
od:
od:
plot3d( blocks, s=0..1, t=0..1, grid=[2,2]) ;  end:

mvcal[crossproduct] :=
  proc(one, two, initpoint)
  local result, point;
  if not type(one, '[algebraic, algebraic, algebraic]') or
     not type(two, '[algebraic, algebraic, algebraic]') 
      then
      ERROR(`Illegal formats`);
  fi;
  result := linalg[crossprod](one, two);
  print(result);  point := op(2, initpoint) ;
  plot3d({[point[1] + one[1]*t, point[2] + one[2]*t, point[3] + one[3]*t],
[point[1] + one[1]- s*1/2*one[1],point[2] + one[2]-s*one[2]/2, point[3] +
one[3]-s*one[3]],
[point[1] + one[1]-s*one[1], point[2] + one[2]-s*one[2], point[3] +
one[3]-s*one[3]/2],
          [point[1] + two[1]*t, point[2] + two[2]*t, point[3] + two[3]*t],
[point[1] + two[1]- s*1/2*two[1], point[2] + two[2]-s*two[2]/2, point[3] +
two[3]-s*two[3]],
[point[1] + two[1]-s*two[1], point[2] + two[2]-s*two[2], point[3] +
two[3]-s*two[3]/2],
          [point[1] + result[1]*t, point[2] + result[2]*t, point[3] +
result[3]*t],
[point[1] + result[1]- s*1/2*result[1], point[2] + result[2]-s*result[2]/2,
point[3] + result[3]-s*result[3]],
[point[1] + result[1]-s*result[1], point[2] + result[2]-s*result[2], point[3] +
result[3]-s*result[3]/2]},
          t = 0..1, s = 0..0.3, axes=NORMAL, grid=[2,2], scaling=CONSTRAINED);
end:

mvcal[cyldemo] := proc(theta, radius, height) ;  
if (radius <0 ) then ERROR(`invalid input`) fi ; 
if (evalf(theta) < 0  or  evalf(theta - 2*Pi) > 0) then ERROR(`invalid input`) 
fi ; 
t := 't'; u:= 'u'; v:= 'v' ;  
plots[animate3d]({[ ((radius-1)*(max(1,min(2, t))-1) +
1)*cos(min(1,t)*theta)*u,  
((radius-1)*(max(1,min(2, t))-1) + 1)*sin(min(1,t)*theta) *u,
height*(max(2,min(3,t))-2)], [((radius-1)*(max(1,min(2, t))-
1) + 1)*cos(min(1,t)*theta) +0.05*radius*u, ((radius-1)*(max(1,min(2, t))-1) +
1)*sin(min(1,t)*theta) +0.05*radius*v, height*(max(2,min(3,t))-2)]},  u=0..1, 
v=0..1, t=0..3, colour= RED, grid=[2, 2], frames = 30, orientation=[30, 70],
axes= NORMAL, scaling = CONSTRAINED) ;  end:

mvcal[drdtplot] := proc (rboundary , trange)
    local line, f, g, a, b, i;
    f := op(1, op(2, rboundary));
    g := op(2, op(2, rboundary));
    a := op(1, op(2, trange));
    b := op(2, op(2, trange));
    theta := 'theta' ;
if not type (trange, string = range) or not (op(1, rboundary) = 'r') or not
(op(1, trange) = 'theta') or not type (evalf(a), numeric) or not type(evalf(b),
numeric) then
ERROR(`input expression is not of drdt type` ) fi; 
    line := {[f*cos(theta), f*sin(theta),theta=a..b],
                 [g*cos(theta), g*sin(theta),theta=a..b],[0,0,0,0]};
    for i from 0 to 5 do
        if evalf( subs(theta=(5-i)/5*a + i/5*b,g)) - evalf(subs(theta=
(5-i)/5*a + i/5*b,f)) < -0.0001
        then ERROR(`the positions of the radiall functions are not correct`)
fi;
    line := line union {[subs(theta=evalf((5-i)/5*a + i/5*b),f*cos(theta)),
 subs(theta=evalf((5-i)/5*a + i/5*b),f*sin(theta)),subs(theta=evalf((5-i)/5*a +
i/5*b),g*cos(theta)),subs(theta=evalf((5-i)/5*a + i/5*b),g*sin(theta))]};
    od;
    plot(line);
    end:

mvcal[dtdrplot] := proc (tbound , rrange)
    local graph1, line , graph2, f,g,a,b, h1, h2, h, i;
        f:= op( 1, op(2,tbound));
        g := op(2, op(2, tbound));
        a := op(1, op(2, rrange));
        b := op(2, op(2, rrange));
        r := 'r' ; s := 's' ; 
if not type (rrange, string = range) or not (op(1, tbound) = 'theta') or not
(op(1, rrange) = 'r') or not type (evalf(a), numeric) or not type(evalf(b),
numeric) then
ERROR(`input expression is not of dtdr type` ) fi; 
    graph1 := plot ({[r*cos(f), r*sin(f), r=a..b],
                 [r*cos(g), r*sin(g), r=a..b], [0,0,0,0]}) :
line := {} :
    for i from 0 to 5 do
		h := (5-i)/5*a + i/5*b ;
        h1 := evalf(subs(r=h,g)) : 
		h2 := evalf( subs( r= h ,f)) :
if (h1 - h2 < - 0.000001)
        then ERROR(`the positions of the angle functions are not correct`)  fi;
		if (h1 -h2 > 10^(-3) ) then 
        line := line union {[h*cos(s), h*sin(s), s= h2..h1]};  fi ; od;
		graph2 := plot(line) :
    plots[display] ( {graph1, graph2 } ) ; 
    end:

mvcal[dxdydzplot] := proc(xbound, ybound, zrange)
    local surface,f,g,h,k,a,b ;
        f := op(1, op(2, xbound));
        g := op(2, op(2, xbound));
        h := op(1, op(2, ybound));
        k := op(2, op(2, ybound));
        a := op(1, op(2, zrange));
        b := op(2, op(2, zrange));
	if not type (zrange, string = range) or not (op(1, xbound) = 'x') or not
(op(1, ybound) = 'y') or not (op(1, zrange) = 'z') or not type (evalf(a),
numeric) or not type(evalf(b), numeric) then
ERROR(`input expression is not of dxdydz type` ) fi; 			z := 'z' :
        y := 'y' ;
    surface := {[subs(z=u*a+(1-u)*b, y=subs(z=u*a+(1-u)*b, v*h+(1-v)*k),f), 
subs(
z=u*a+(1-u)*b,v*h+(1-v)*k), u*a+(1-u)*b],
[ subs(z=u*a+(1-u)*b, y=subs(z=u*a+(1-u)*b, v*h+(1-v)*k),g),
subs(z=u*a+(1-u)*b,
v*h+(1-v)*k) , u*a+(1-u)*b],
[v*subs(z=u*a+(1-u)*b,y=subs(z=u*a+(1-u)*b,
 h),f) + (1-v)*subs(z=u*a+(1-u)*b,y=subs(z=u*a+(1-u)*b, h), g),
subs(z=u*a+(1-u)
*b,h), u*a+(1-u)*b],
[  v*subs(z=u*a+(1-u)*b,y=subs(z=u*a+(1-u)*b,
 k),f) + (1-v)*subs(z=u*a+(1-u)*b,y=subs(z=u*a+(1-u)*b, k), g),
subs(z=u*a+(1-u)
*b,k), u*a+(1-u)*b],
[  v*subs(z=a, y=u*subs(z=a,h)+ (1-u)*subs(z=a,k), f) + (1-v)*subs(z=a,
y=u*subs(z=a,h)+ (1-u)*subs(z=a,k), g), u*subs(z=a,h)+ (1-u)*subs(z=a,k), a],
[v*subs(z=b, y=u*subs(z=b,h)+ (1-u)*subs(z=b,k), f) + (1-v)*subs(z=b,
y=u*subs(z=b,h)+ (1-u)*subs(z=b,k), g), u*subs(z=b,h)+ (1-u)*subs(z=b,k) , b]};
    plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[15,15]);
end :

mvcal[dxdyplot] := proc(xboundary, yrange)
    local line ,i,f,g,a,b;
    f:= op(1,op(2,xboundary));
    g:= op(2,op(2,xboundary));
    a:= op(1,op(2,yrange));
    b:= op(2, op(2, yrange));
    y:= 'y' ;
if not type (yrange, string = range) or not (op(1, xboundary) = 'x') or not
(op(1, yrange) = 'y') or not type (evalf(a), numeric) or not type(evalf(b),
numeric) then
ERROR(`input expression is not of dxdy type` ) fi; 
    line := {[f,y,yrange],[g,y,yrange] , [0,0,0,0]};
    for i from 0 to 5 do
        if evalf(subs(y=(5-i)/5*a+i/5*b,g)) - evalf( subs(y=(5-i)/5*a+i/5*b,f))
< -0.000001
        then ERROR(`the positions of the functions are not correct`) fi;
        line := line union {[subs(y=evalf((5-i)/5*a+i/5*b),f),
evalf((5-i)/5*a+i
/5*b),
    subs(y=evalf((5-i)/5*a+i/5*b),g),evalf((5-i)/5*a+i/5*b)]};
        od;
plot(line);
end:

mvcal[dxdzdyplot] := proc(xbound, zbound, yrange)
    local surface,f,g,h,k,a,b;
        f:= op(1, op(2, xbound));
        g:= op(2, op(2, xbound));
        h:= op(1, op(2, zbound));
        k:= op(2, op(2, zbound));
        a:= op(1, op(2, yrange));
        b:= op(2, op(2, yrange));
	if not type (yrange, string = range) or not (op(1, xbound) = 'x') or not
(op(1, zbound) = 'z') or not (op(1, yrange) = 'y') or not type (evalf(a),
numeric) or not type(evalf(b), numeric) then
ERROR(`input expression is not of dxdzdy type` ) fi; 			y := 'y' :
         z := 'z' :
    surface := {[subs(y=u*a+(1-u)*b, z=subs(y=u*a+(1-u)*b,
v*h+(1-v)*k),f),u*a+(
1-u)*b, subs(y=u*a+(1-u)*b,v*h+(1-v)*k)],
[ subs(y=u*a+(1-u)*b, z=subs(y=u*a+(1-u)*b, v*h+(1-v)*k),g), u*a+(1-u)*b,
subs(y
=u*a+(1-u)*b,v*h+(1-v)*k) ],
[v*subs(y=u*a+(1-u)*b,z=subs(y=u*a+(1-u)*b,
 h),f) + (1-v)*subs(y=u*a+(1-u)*b,z=subs(y=u*a+(1-u)*b, h), g), u*a+(1-u)*b,
subs(
y=u*a+(1-u)*b,h)],
[  v*subs(y=u*a+(1-u)*b,z=subs(y=u*a+(1-u)*b,
 k),f) + (1-v)*subs(y=u*a+(1-u)*b,z=subs(y=u*a+(1-u)*b, k), g),
u*a+(1-u)*b,subs
(y=u*a+(1-u)*b,k)],
[  v*subs(y=a, z=u*subs(y=a,h)+ (1-u)*subs(y=a,k), f) + (1-v)*subs(y=a,
z=u*subs(y=a,h)+ (1-u)*subs(y=a,k), g), a, u*subs(y=a,h)+ (1-u)*subs(y=a,k)],
[v*subs(y=b, z=u*subs(y=b,h)+ (1-u)*subs(y=b,k), f) + (1-v)*subs(y=b,
z=u*subs(y=b,h)+ (1-u)*subs(y=b,k), g), b, u*subs(y=b,h)+ (1-u)*subs(y=b,k) ]};
    plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[15,15]);
end:

mvcal[dydxdzplot] := proc(ybound, xbound, zrange)
    local surface,f,g,h,k,a,b;
        f:= op(1, op(2, ybound));
        g:= op(2, op(2, ybound));
        h:= op(1, op(2, xbound));
        k:= op(2, op(2, xbound));
        a:= op(1, op(2, zrange));
        b:= op(2, op(2, zrange));
if not type (zrange, string = range) or not (op(1, xbound) = 'x') or not (op(1,
ybound) = 'y') or not (op(1, zrange) = 'z') or not type (evalf(a), numeric) or
not type(evalf(b), numeric) then
ERROR(`input expression is not of dydxdz type` ) fi; 
        x := 'x' :
        z := 'z' :
    surface := {[subs(z=u*a+(1-u)*b,v*h+(1-v)*k),
subs(z=u*a+(1-u)*b, x=subs(z=u*a+(1-u)*b, v*h+(1-v)*k),f),u*a+(1-u)*b],
[ subs(z=u*a+(1-u)*b,v*h+(1-v)*k), subs(z=u*a+(1-u)*b,
x=subs(z=u*a+(1-u)*b, v*h+(1-v)*k),g), u*a+(1-u)*b],
[subs(z=u*a+(1-u)*b,h), v*subs(z=u*a+(1-u)*b,x=subs(z=u*a+(1-u)*b,
 h),f) + (1-v)*subs(z=u*a+(1-u)*b,x=subs(z=u*a+(1-u)*b, h), g), u*a+(1-u)*b],
[  subs(z=u*a+(1-u)*b,k),v*subs(z=u*a+(1-u)*b,x=subs(z=u*a+(1-u)*b,
 k),f) + (1-v)*subs(z=u*a+(1-u)*b,x=subs(z=u*a+(1-u)*b, k), g), u*a+(1-u)*b],
[  u*subs(z=a,h)+ (1-u)*subs(z=a,k), v*subs(z=a,
x=u*subs(z=a,h)+ (1-u)*subs(z=a,k), f) + (1-v)*subs(z=a,
x=u*subs(z=a,h)+ (1-u)*subs(z=a,k), g), a],
[u*subs(z=b,h)+ (1-u)*subs(z=b,k), v*subs(z=b,
x=u*subs(z=b,h)+ (1-u)*subs(z=b,k), f) + (1-v)*subs(z=b,
x=u*subs(z=b,h)+ (1-u)*subs(z=b,k), g), b]};
    plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[15,15]);
end:

mvcal[dydxplot] := proc(ybound , xrange)
    local line, f,g,a,b ,i;
    f := op(1, op(2,ybound));
    g := op(2, op(2,ybound));
    a := op(1, op(2,xrange));
    b := op(2, op(2,xrange));
    x := 'x';
if not type(xrange, algebraic = range) or not (op(1,ybound) = 'y') or not
(op(1, xrange) = 'x' ) or not type(evalf(a ), numeric) or not type(evalf(b) ,
numeric) then ERROR(`Input expression is not of dydx type`)  fi;
 line := {[x,f,x=a..b],[x,g,x=a..b],[0,0,0,0]};
    for i from 0 to 5 do
        if ( evalf(subs(x=(5-i)/5*a+i/5*b,g)) - evalf(
subs(x=(5-i)/5*a+i/5*b,f)
) < -0.00001)  then ERROR(`arrangements of the functions are incorrect`)
        fi;
      line := line union
{[evalf((5-i)/5*a+i/5*b),subs(x=evalf((5-i)/5*a+i/5*b),
 f), evalf((5-i)/5*a+i/5*b), subs(x=evalf((5-i)/5*a+i/5*b),g)]};
        od;
plot(line);
end:

mvcal[dydzdxplot] := proc(ybound, zbound, xrange)
    local surface,f,g,h,k,a,b ;
        f:= op(1, op(2, ybound));
        g:= op(2, op(2, ybound));
        h:= op(1, op(2, zbound));
        k:= op(2, op(2, zbound));
        a:= op(1, op(2, xrange));
        b:= op(2, op(2, xrange));
if not type (xrange, string = range) or not (op(1, ybound) = 'y') or not (op(1,
zbound) = 'z') or not (op(1, xrange) = 'x') or not type (evalf(a), numeric) or
not type(evalf(b), numeric) then
ERROR(`input expression is not of dydzdx type` ) fi; 
         x := 'x' :
         z := 'z' :
    surface := {[u*a+(1-u)*b, subs(x=u*a+(1-u)*b, z=subs(x=u*a+(1-u)*b,
v*h+(1-v
)*k),f), subs(x=u*a+(1-u)*b,v*h+(1-v)*k)],
[u*a+(1-u)*b, subs(x=u*a+(1-u)*b, z=subs(x=u*a+(1-u)*b, v*h+(1-v)*k),g), 
subs(x
=u*a+(1-u)*b,v*h+(1-v)*k) ],
[u*a+(1-u)*b, v*subs(x=u*a+(1-u)*b,z=subs(x=u*a+(1-u)*b,
 h),f) + (1-v)*subs(x=u*a+(1-u)*b,z=subs(x=u*a+(1-u)*b, h), g), 
subs(x=u*a+(1-u
)*b,h)],
[ u*a+(1-u)*b, v*subs(x=u*a+(1-u)*b,z=subs(x=u*a+(1-u)*b,
 k),f) + (1-v)*subs(x=u*a+(1-u)*b,z=subs(x=u*a+(1-u)*b, k), g),
subs(x=u*a+(1-u)
*b,k)],
[  a,  v*subs(x=a, z=u*subs(x=a,h)+ (1-u)*subs(x=a,k), f) + (1-v)*subs(x=a,
z=u*subs(x=a,h)+ (1-u)*subs(x=a,k), g),u*subs(x=a,h)+ (1-u)*subs(x=a,k)],
[b, v*subs(x=b, z=u*subs(x=b,h)+ (1-u)*subs(x=b,k), f) + (1-v)*subs(x=b,
z=u*subs(x=b,h)+ (1-u)*subs(x=b,k), g),  u*subs(x=b,h)+ (1-u)*subs(x=b,k) ]};
    plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[15,15]);
end:

mvcal[dzdxdyplot] := proc(zbound, xbound, yrange)
    local surface,f,g,h,k,a,b ;
        f:= op(1, op(2, zbound));
        g:= op(2, op(2,  zbound));
        h:= op(1, op(2, xbound));
        k:= op(2, op(2, xbound));
        a:= op(1, op(2, yrange));
        b:= op(2, op(2, yrange));
if not type (yrange, string = range) or not (op(1, xbound) = 'x') or not (op(1,
zbound) = 'z') or not (op(1, yrange) = 'y') or not type (evalf(a), numeric) or
not type(evalf(b), numeric) then
ERROR(`input expression is not of dzdxdy type` ) fi; 
        y := 'y' :
        x := 'x' :
   surface := {[subs(y=u*a+(1-u)*b,v*h+(1-v)*k),u*a+(1-u)*b,
subs(y=u*a+(1-u)*b, x=subs(y=u*a+(1-u)*b, v*h+(1-v)*k),f)],
[ subs(y=u*a+(1-u)*b,v*h+(1-v)*k), u*a+(1-u)*b, subs(y=u*a+(1-u)*b,
x=subs(y=u*a+(1-u)*b, v*h+(1-v)*k),g)],
[subs(y=u*a+(1-u)*b,h), u*a+(1-u)*b, v*subs(y=u*a+(1-u)*b,x=subs(y=u*a+(1-u)*b,
 h),f) + (1-v)*subs(y=u*a+(1-u)*b,x=subs(y=u*a+(1-u)*b, h), g)],
[  subs(y=u*a+(1-u)*b,k),
u*a+(1-u)*b,v*subs(y=u*a+(1-u)*b,x=subs(y=u*a+(1-u)*b,
 k),f) + (1-v)*subs(y=u*a+(1-u)*b,x=subs(y=u*a+(1-u)*b, k), g)],
[  u*subs(y=a,h)+ (1-u)*subs(y=a,k), a, v*subs(y=a,
x=u*subs(y=a,h)+ (1-u)*subs(y=a,k), f) + (1-v)*subs(y=a,
x=u*subs(y=a,h)+ (1-u)*subs(y=a,k), g)],
[u*subs(y=b,h)+ (1-u)*subs(y=b,k), b, v*subs(y=b,
x=u*subs(y=b,h)+ (1-u)*subs(y=b,k), f) + (1-v)*subs(y=b,
x=u*subs(y=b,h)+ (1-u)*subs(y=b,k), g)]};
    plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[15,15]);
end:

mvcal[dzdydxplot] := proc(zbound,ybound,xrange)
    local surface,f,g,h,k,a,b;
        f:= op(1, op(2, zbound));
        g:= op(2, op(2, zbound));
        h:= op(1, op(2, ybound));
        k:= op(2, op(2, ybound));
        a:= op(1, op(2, xrange));
        b:= op(2, op(2, xrange));
if not type (xrange, string = range) or not (op(1, ybound) = 'y') or not (op(1,
zbound) = 'z') or not (op(1, xrange) = 'x') or not type (evalf(a), numeric) or
not type(evalf(b), numeric) then
ERROR(`input expression is not of dzdydx type` ) fi; 
        y := 'y' :
        x := 'x' :
    surface := {[u*a+(1-u)*b,subs(x=u*a+(1-u)*b,v*h+(1-v)*k),
subs(x=u*a+(1-u)*b, y=subs(x=u*a+(1-u)*b, v*h+(1-v)*k),f)],
[ u*a+(1-u)*b, subs(x=u*a+(1-u)*b,v*h+(1-v)*k), subs(x=u*a+(1-u)*b,
y=subs(x=u*a+(1-u)*b, v*h+(1-v)*k),g)],
[u*a+(1-u)*b,subs(x=u*a+(1-u)*b,h), v*subs(x=u*a+(1-u)*b,y=subs(x=u*a+(1-u)*b,
 h),f) + (1-v)*subs(x=u*a+(1-u)*b,y=subs(x=u*a+(1-u)*b, h), g)],
[u*a+(1-u)*b,subs(x=u*a+(1-u)*b,k),v*subs(x=u*a+(1-u)*b,y=subs(x=u*a+(1-u)*b,
 k),f) + (1-v)*subs(x=u*a+(1-u)*b,y=subs(x=u*a+(1-u)*b, k), g)],
[ a, u*subs(x=a,h)+ (1-u)*subs(x=a,k), v*subs(x=a,
y=u*subs(x=a,h)+ (1-u)*subs(x=a,k), f) + (1-v)*subs(x=a,
y=u*subs(x=a,h)+ (1-u)*subs(x=a,k), g)],
[b,u*subs(x=b,h)+ (1-u)*subs(x=b,k), v*subs(x=b,
y=u*subs(x=b,h)+ (1-u)*subs(x=b,k), f) + (1-v)*subs(x=b,
y=u*subs(x=b,h)+ (1-u)*subs(x=b,k), g)]};
    plot3d(surface, u=0..1, v=0..1, axes=BOXED, grid =[15,15]);
end:

mvcal[lagrange] := proc(f,c,g,h,v)
        local curve;
                curve := solve(f=c, y);
        plot({curve, solve(g,y)},h,v);
end:

mvcal[levelcurve] := proc(F, c, xrange, yrange)
local a1, a2,a3, a4, b1, b2, b3,b4,c1, e,f, curve,r, i ,ans, j;
        if type(c, constant) then c1 := [c] else c1 := c fi;
        if not type(c1, list(constant)) then ERROR(`illegal constant`) fi;
        curve := {};
        b3 := diff(diff(diff(F,y),y),y) ;
        a3 := diff(diff(diff(F,x),x),x) ;
        if b3 =0 then
                if a3 =0 then
                        for i from 1 to nops(c1) do
                                ans := {solve(F=c1[i],x)};
                                for j from 1 to nops(ans) do
                                        curve := curve union { [subs(y=t,
ans[j]), t, t =op(2,yrange)]};
                                od;
                        od;
                        for i from 1 to nops(c1) do
                                ans := {solve(F=c1[i],y)};
                                for j from 1 to nops(ans) do
                                        curve := curve union { [t, subs(x=t,
ans[j]), t=op(2,xrange)]};
                                od;
                        od;
                        plot(curve, xrange, yrange);
                else
                        for i from 1 to nops(c1) do
                        curve := curve union {solve(F=c1[i],y)};
                        od;
                plot(curve, xrange, yrange);
                fi;
        else
                if a3 =0 then
                        for i from 1 to nops(c1) do
                                ans := {solve(F=c1[i],x)};
                                for j from 1 to nops(ans) do
                                        curve := curve union { [ans[j], y,
yrange]};
                                        od;
                                od;
                                plot(curve, xrange, yrange);
                else
                        if diff( subs(x^4=r, subs(x^2=r, f)),x) = 0 then
                                for i from 1 to nops(c1) do
                                ans := {solve(F=c1[i],x)};
                                for j from 1 to nops(ans) do
                                        curve := curve union { [ans[j], y,
yrange]};
                                        od;
                                od;
                                plot(curve, xrange, yrange);
                        else
                                for i from 1 to nops(c1) do
                                curve := curve union {solve(F=c1[i],y)};
                                od;
                        plot(curve, xrange, yrange);
                        fi;
                fi;
        fi;
        end :

mvcal[levelsurface] := proc(f, c, h, v)
        local surface, i;
       if nops(c) = 1 then
                surface := solve(f = c, z);
                plot3d({surface}, h, v);
        else
                surface := {};
                for i from 1 to nops(c) do
                        surface := surface union {solve(f = c[i], z)};
                od;
                plot3d(surface, h, v);
        fi;
end :

mvcal[plotsvf3d]  := proc (G, F, urange, vrange)
local a,b,c,d,u1,v1, Gsubs, Fsubs,ud,vd,points,graph1,graph2;
if type(urange, algebraic=range) and type(vrange, algebraic=range) then
    a := evalf(op(1, op(2, urange)));    b := evalf(op(2, op(2, urange)));
    c := evalf(op(1, op(2, vrange)));    d := evalf(op(2, op(2, vrange)));
    if ((a>b) or (c>d)) then ERROR(`Bad range value`) fi;
   else ERROR(`Illegal range format`) fi;
ud := (b-a)/6;   vd := (d-c)/6;
points :={} ;
    for u1 from a by ud to b  do
        for v1 from c by vd to d do
			Gsubs := subs(u=u1, v=v1, G) :  
			Fsubs := subs(x=Gsubs[1], y=Gsubs[2], z=Gsubs[3], F) :
			points := points union {[Gsubs[1]+t*Fsubs[1], Gsubs[2]+t*Fsubs[2],
Gsubs[3]+t*Fsubs[3]] } ;  od ; od;
with (plots):
graph1 := plot3d(points, t=0..1, s=0..1, grid=[2,2]):
graph2 := plot3d( G, urange, vrange); 
display3d({graph1, graph2});
end:


mvcal[plotvectorfield]  := proc(F, xrange, yrange, increment)
    local  func,t,xlabel,ylabel,
           curvseq,xmin,xmax,ymin,ymax,xx,yy,vector,m,magn,nx,ny,vlen,label;
   if type(xrange, algebraic = range) then
        xlabel := LABEL(TEXT(convert(op(1,xrange),name)));
        xmin := evalf(op(1,op(2,xrange)));
        xmax := evalf(op(2,op(2,xrange)));
        if not member(x, {op(1,xrange)}) then
            print(`warning:  x not used in x range`)
        fi;
    else ERROR(`invalid x range`)
    fi;
    if type(args[3], algebraic = range) then
        ylabel := LABEL(TEXT(convert(op(1,args[3]),name)));
        ymin := evalf(op(1,op(2,args[3])));
        ymax := evalf(op(2,op(2,args[3])));
        if not member(y, {op(1,yrange)}) then
            print(`warning:  y not used in y range`)
        fi;
    else ERROR(`invalid y range`)
    fi;
   curvseq := NULL;
    fb := 1;
   func := F;
     x1 := xmin; x2 := xmax;
    y1 := ymin; y2 := ymax;
	curvseq := { } : 
    magn := (xmax - xmin)/50;
vector1 := subs(x=x1, y=y1, F) :  vector2 := subs(x=x1,y=y2, F) :  vector3 :=
subs(x=x2,y=y1, F) : vector4 := subs(x=x2, y=y2, F) ;
vlen := sqrt((vector1[1])^2 +(vector1[2])^2+ (vector2[1])^2 +(vector2[2])^2
+(vector3[1])^2 +(vector3[2])^2+(vector4[1])^2+(vector4[2])^2)/4 ;
if vlen = 0 then
                m[1] := 0;  m[2] := 0
            else
                m[1] := magn/vlen;  m[2] := magn/vlen
            fi;
     for xx from xmin by evalf(increment) to xmax do
        for yy from ymin by evalf(increment) to ymax do
          nx := subs(x=xx, y=yy, F)[1] ;  ny := subs(x=xx, y=yy, F)[2] ;       
   
            curvseq := curvseq union { [xx,yy,xx+nx,yy+ny],
[xx+nx-m[1]*nx-m[2]*ny,yy+ny-m[2] *ny+m[1]*nx,xx+nx, 
yy+ny,xx+nx-m[1]*nx+m[2]*ny ,yy+ny-m[2]*ny-m[1]*nx] };
        od;
    od;
x1 := x1 + min(0, vector1[1], vector2[1]) ;  x2 := x2 + max(0, vector3[1],
vector4[1]);  y1 := y1 + min(0, vector1[2], vector3[2]);  y2 := y2 +max(0,
vector2[2], vector4[2]) ;
    plot(curvseq,  x=x1..x2, y=y1..y2, numpoints = 2);
end :

mvcal[plotvectorfield3d] :=
proc (F, xrange, yrange, increment, zrange, zlevel)
local a,b,c,d,e,f,i, j, k, xd,yd,zd,label,xscale, yscale, zscale, zlayers,
test;
if type(xrange, algebraic=range) and type(yrange, algebraic=range)
  and type(zrange, algebraic=range) then
    a := op(1, op(2, xrange));    b := op(2, op(2, xrange));
    c := op(1, op(2, yrange));    d := op(2, op(2, yrange));
    e := op(1, op(2, zrange));    f := op(2, op(2, zrange));
    if (a>b) or (c>d) or (e>f) then ERROR(`Bad range value`) fi;
    label := [convert(op(1,xrange),name),
              convert(op(1,yrange),name),
              convert(op(1,zrange),name)];
  else ERROR(`Illegal range format`)
fi;

if type(zlevel, string=numeric) then
    test := op(1, zlevel);
    if not (test = 'layers') then ERROR(`Unrecognized \'layers\'`)
      else zlayers := eval(op(2, zlevel)) fi;
else ERROR(`Unrecognized layer parameter`);
fi;
xd := eval(((b-a)/increment));  b := a + xd*increment;
yd := eval(((d-c)/increment));  d := c + yd*increment;
zd := eval(zlayers-1);
xscale := (b-a)/xd;    yscale := (d-c)/yd;    zscale := (f-e)/zd;
points := {} ; 
for k from e by zscale to f do
	points := points union { [t*a+(1-t)*b, s*c+(1-s)*d, k]};
	for i from a by increment to b do
		for j from c by increment to d do
		
		u1 := subs(x=i, y=j, z=k, op(1, F)) ; 
		u2 := subs(x=i, y=j, z=k,  op(2, F));  
		u3 := subs(x=i, y=j, z=k,  op(3, F));
		points := points union {[i+t*u1, j+t*u2, k+t*u3]  }; od ; od ; od; 


plot3d(points, t = 0..1, s=0..1, grid=[2,2], orientation=[30,70], style =
WIREFRAME);
end :

mvcal[sphdemo] := proc( phi, theta, rho) ;  
u := 'u'; v:= 'v';
 if (evalf(rho) < 0 ) then ERROR (`invalid input`) fi;
 if evalf(phi) <0 or evalf(phi-Pi) > 0 then ERROR (`invalid input`) fi;
  plots[animate3d]({ [ u*( (rho-1)*(max(2, min(3, t)) - 2) +1)*sin( min(1,
t)*phi) *
 cos( (max( 0, min(2, t) -1))*theta),
 u*( (rho-1)*(max(2, min(3, t)) - 2) +1)*sin( min(1, t)*phi) * sin( (max( 0,
 min(2, t) -1))*theta) ,
  u* ((rho-1)*(max(2, min(3, t)) - 2) +1)*cos( min(1, t)*phi) ],
   [ ( (rho-1)*(max(2, min(3, t)) - 2) +1)*sin( min(1, t)*phi) * cos( (max( 0,
 min(2, t) -1))*theta) +rho*0.05*v,
 ((rho-1)*(max(2, min(3, t)) - 2) +1)*sin( min(1, t)*phi) * sin( (max( 0,
min(2,t) -1))*theta) +rho*0.05*u ,
 ((rho-1)*(max(2, min(3, t)) - 2) +1)*cos( min(1, t)*phi) ]},
 u=0..1, v=0..1, t=0..3,  frames= 30,  grid=[2,2], scaling = CONSTRAINED,
 axes = NORMAL, colour= RED, orientation = [30, 70]) ; end :

mvcal[ysliceapp] := proc (f, xexpress, yexpress) 
local slices, x1, x2, a, b, c, d ;
        a := op(1, op(2, xexpress));
        b := op(2, op(2, xexpress));
        c := op(1, op(2, yexpress));
        d := op(2, op(2, yexpress));
x := 'x' : y := 'y' :
slices := {[b,c*t+(1-t)*d,s*(subs(x=b,y=c*t+(1-t)*d, f))]}:
for i from 1 to 5 do
x1 := a+(i-1)*(b-a)/5 :
x2 := a+ i*(b-a)/5 :
slices := slices union {[s*x1+(1-s)*x2, c*t+(1-t)*d ,subs(x= x2, y
=c*t+(1-t)*d, f)],
[s*x1+(1-s)*x2, c,  t* subs(x=x2 , y=c, f)], [s*x1+(1-s)*x2, d, t*subs(x=x2 ,
y=d, f)]} :
od :
plot3d (slices, s=0..1, t=0..1, grid=[2,20]);
end :

mvcal[demo] := proc (numbers) ;                      
if (numbers = 12.1) then plot3d({[t,t,t], [1/2*sin(t*Pi/4), 1/2*sin(t*Pi/4),
1/2*cos(t*Pi/4)]}, t=0..1, s
=0..1,
grid =[35,2], orientation=[20,80], axes=NORMAL) ; 

elif (numbers = 27.1) then plot3d({4*y^2, 0}, x=0..2, y=-1..1);

elif(numbers=27.2) then 
plot3d({[r*cos(theta), r*sin(theta), 4 -2*r],[r*cos(theta), 
r*sin(theta),0]}, r=1..1.5, theta=0..2*Pi);

elif(numbers=27.3) then
 plot3d( {[t, s*(3-t)+(1-s)*(2*t),0],[t, s*(3-t)+(1-s)*(2*t), 
1/2(s*(3-t)+(1-s)*(2*t))^2+ t*s*(3-t)+(1-s)*(2*t)+ 3*cos(6*t)]} , 
t=0..1, s=0..1) ; 

elif (numbers = 15.1) then
f := (y^3 - 3*y)/(1+x^2) :
plot3d( f, x=-2..2, y=-2..2)  ;

elif( numbers = 15.2) then
fig2 := plot3d( { [t,1/2,subs(x=t,y=1/2,f)], [0, (t+1/2), subs(x=0, y=t+1/2,
f)], [t,t+1/2,subs(x=t,y=t+1/2,f)],[ -t, t+1/2, subs(x=-t,y=t+1/2,f)],
[2*t,t+1/2,subs(x=2*t, y=t+1/2,f)], [-1/2*t, t+1/2, subs(x=-1/2*t, y=t+1/2,
f)]}, t=-1..1, s=0..1, grid=[30,2], scaling = CONSTRAINED) :
normal :=  plot3d ({[0, 1/2 + 9/4*t, -11/8 + t], [0,11/4-0.5*t, -3/8-
0.4*t],[0,11/4-0.5*t,-3/8-0.1*t]}, t=0..1, s=0..2, grid=[2,2]) :
with(plots):
display3d({fig2}) ;

elif (numbers=15.3) then

display3d({fig2, normal}) ;  

elif (numbers = 15.4) then
graph := plot3d(f,x=-1..1, y=-0.5..1.1) :
tangent := plot3d( -9/4*y + 9/4/2-11/8, x=-1..1, y=0..1, scaling =
CONSTRAINED):
normal :=  plot3d ({[0, 1/2 + 9/4*t, -11/8 + t], [0,11/4-0.5*t,
-3/8-0.4*t],[0,11/4-0.5*t,-3/8-0.1*t]}, t=0..1, s=0..2, grid=[2,2]) :
with(plots): 
display3d( {normal, tangent});    

elif(numbers = 15.5) then
display3d({graph, tangent}) ;

elif(numbers = 15.6) then
display3d( {graph, tangent, normal}) ;

elif (numbers = 21.1) then
graph1 := plot3d( x^2 + y^2+2, x=-1..1, y= -1..1 ):
graph2 := plot3d( 10 - (x-1)^2 - y^2, x=-1..1, y= -1..1):
graph3 := plot3d([1/3*(2-cos(2*(t+2*Pi/3)))*cos(t),
1/4*(2-cos(2*t))*sin(t),s], t=0..2*Pi, s=0..11): 
plots[display3d]({graph1, graph2, graph3}); 

elif (numbers = 21.2)  then
graph1 := plot3d( [x^2 + y^2+2+5,x,y], x=-1..1, y= -1..1 ):
graph2 := plot3d( [18 - (x-1)^2 - y^2, x,y], x=-1..1, y= -1..1):
graph3 := plot3d([s, 1/3*(2-cos(2*(t+2*Pi/3)))*cos(t),
1/4*(2-cos(2*t))*sin(t)], t=0..2*Pi, s=5..20):
graph4 := plot3d([0, 1/3*(2-cos(2*(t+2*Pi/3)))*cos(t),
1/4*(2-cos(2*t))*sin(t)], t=0..2*Pi, s=0..1, grid=[25,2]):
plots[display3d]({graph1, graph2, graph3, graph4}); 

elif (numbers = 21.3)   then
graph1 := plot3d( [x, x^2 + y^2+2+5,y], x=-1..1, y= -1..1 ):
graph2 := plot3d( [x, 18 - (x-1)^2 - y^2,y], x=-1..1, y= -1..1):
graph3 := plot3d([ 1/3*(2-cos(2*(t+2*Pi/3)))*cos(t),s,
1/4*(2-cos(2*t))*sin(t)], t=0..2*Pi, s=5..20):
graph4 := plot3d([ 1/3*(2-cos(2*(t+2*Pi/3)))*cos(t),0,
1/4*(2-cos(2*t))*sin(t)], t=0..2*Pi, s=0..1, grid=[25,2]):
plots[display3d]({graph1, graph2, graph3, graph4});  

elif (numbers = 21.4)   then
graph1 := plot3d(x^2+y^2, x=-2..2,y=-2..2) :
graph2 := plot3d(4-x^2 -y^2, x=-2..2, y=-2..2):
plots[display3d]( {graph1, graph2}); 

elif (numbers = 21.5)  then
graph1 := plot3d(x^2+y^2, x=-2..2,y=-2..2) :
graph2 := plot3d(4-x^2 -y^2, x=-2..2, y=-2..2):
graph3 := plot3d([sqrt(2)*cos(theta), sqrt(2)*sin(theta), z], theta=0..2*Pi,
z=-7..7) :
with(plots):
display3d( {graph1, graph2,graph3});    

elif (numbers = 27.1)  then
plot3d({4*y^2, 0}, x=0..2, y=-1..1);
plot3d({[r*cos(theta), r*sin(theta), 4 -2*r],[r*cos(theta),
r*sin(theta),0]}, r=1..1.5, theta=0..2*Pi);
 plot3d( {[t, s*(3-t)+(1-s)*(2*t),0],[t, s*(3-t)+(1-s)*(2*t),
1/2(s*(3-t)+(1-s)*(2*t))^2+ t*s*(3-t)+(1-s)*(2*t)+ 3*cos(6*t)]} ,
t=0..1, s=0..1) ; 

elif (numbers = 3.5)  then
with( mvcal): plotvectorfield ( [1, 1.5], x = -2..2, y=-2..2, 1); 

elif (numbers = 3.6)  then
plot3d( {[(1-t)*2,(1-t)*2,(1-t)*5],[(1-t)*(-1),(1-t)*1,(1-t)*1],
 [(1-t)*1,(1-t)*1,0],[(1-t)*(-1),2*t+(1-t)*1,t+(1-t)*1],[(1-t)*2,2*t+(1-t)*4,
 t+(1-t)*6],[(1-t)*1,2*t+(1-t)*1,t], [3*t+(1-t)*1,3*t+(1-t)*1, 5*t],
 [3*t+(1-t)*2,3*t+(1-t)*4, 5*t+(1-t)*6],[3*t+(1-t)*2,3*t+(1-t)*2, 5*t+(1-t)*5],
 [1*t+(1-t)*2, 3*t+(1-t)*2, 6*t+(1-t)*5], [1*t+(1-t)*2, 3*t+(1-t)*4,
 6*t+(1-t)*6],[1*t+(1-t)*(-1), 3*t+(1-t)*1, 6*t+(1-t)*1]}, t=0..1, s=0..1, grid
=[2,2],
 orientation=[23,73]); 

elif (numbers = 3.1)  then
plot3d([0,0,0],x =-2..2, y=-2..2, grid=[2,2], axes =NORMAL);  

elif (numbers = 3.2)  then
plot3d({[x,0,0], [0,y,0], [0,0,8]}, x=0..4, y=0..6 , grid=[2,2], axes = NORMAL
) ;  

elif (numbers = 3.3)   then
plot3d( [sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)], phi=0..Pi,
theta=0
..2*Pi , scaling  = CONSTRAINED) ;  

elif (numbers = 3.4)   then
plot3d( [sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)], phi=0..Pi,
theta=0
..2*Pi , scaling  = CONSTRAINED) ;  

elif (numbers = 4.1)  then
lines := plot3d( { [0.9*1,0.9*2, 0.9*3], [2*1, 2*2, 2*3], [2.5*1, 2.5*2,
2.5*3],
 [3.2*1, 3.2*2, 3.2*3],
[-1.3*1, -1.3*2, -1.3*3], [-2.2*1, -2.2*2, -2.2*3],[-1.5*1, -1.5*2,
-1.5*3],[-0.5*1, -0.5*2,-0.5*3], [-1.8*1, -1.8*2,-1.8*3],[1.5*1,
1.5*2,1.5*3]
 }, s=0..1, t=0..1, grid=[2,2], style=POINT ) :
vector := plot3d( [t, 2*t, 3*t], t=0..1, s=0..1, style =WIREFRAME, grid=[50,2],
orientation=[15,70]):
with(plots) :
display3d({lines, vector});  

elif (numbers = 4.2)  then
plot3d([t, t*2, t*3], t=-2..2, s=0..1, grid = [2,2], orientation=[15,70]); 

elif (numbers = 4.3)  then
plot3d({[t,t,3*t], [-t,t,0], [0, -3*t, t], [-2*t,-t,t], [5*t,-2*t,-t],
[5*t,4*t,
 -3*t], [-4*t,-5*t, 3*t],
[t, 5*t, -2*t], [0.5*t, -2*t, 0.5*t], [-5*t, 0.5*t, 1.5*t], [-0.5*t, -t,
0.5*t],[t, 2*t,-t], [2*t, -5*t,-1*t], [1,1,3-0.3*t],[1-.2*t, 1 - .2*t,
3-.2*t]},
 t = 0..1,  s=0..1, orientation = [-30,80] , grid=[2,2], scaling =
CONSTRAINED);  else 

fi;
end:


mvcal[ex] := proc( numbers)  
local f, v, line, line1, line2, line3, graph1, tangent, normal, picture,cone,
sphere, graph2, curve1 ; 

if(numbers =4.6) then
plot3d( {[ t+(1-t)*(-2), 2*t+(1-t)*3, 1*t+2*(1-t)],[0, 0,0]}, t=-2..2, s=0..1,
grid=[2,2], axes = NORMAL) ;  

elif (number =4.7) then
plot3d( x+y -1, x=-2..2, y=-2..2,axes=BOXED) ; 

elif (numbers = 4.8)   then
line1 := plot3d( [ 2 -3*t, 3+ 2*t, 1 +t], t=-2..2, s=0..1, grid=[2,2]) :
line2 := plot3d( [2- 2*t, 3+ t, 1-t], t =-2..2, s=0..1, grid =[2,2]) :
line3 := plot3d( [4 + 5*t, 2 -3*t, 2], t=-2..2, s=0..1, grid=[2,2]) :
plots[display3d]({line1, line2, line3}) ;  

elif (numbers = 5.5)  then
graph1 := plot3d( [ 2*sin(t)*cos(t), t, 2*sin(t)*sin(t) ], t=0..2*Pi, s=0..1,
grid = [60,2]) :
graph2 := plot3d( [-sqrt(3)/2 -t, 3*Pi/4 + t, 3/2- sqrt(3)*t], t=0..sqrt(3)/2,
s=0..1, grid=[2,2], axes = NORMAL) : 
plots[display3d] ( {graph1, graph2}, orientation = [51,72]) ;  

elif (numbers = 6.1) then
plot([t, t^2+t+1, t=0..1]);  

elif (numbers = 6.21) then
plot( [3+cos(t), -1+2*sin(t), t=0..2*Pi]) ; 

elif (numbers =6.22) then
plot( [(1+cos(t))*3, (2+sin(t))*2, t=0..2*Pi]) ;

elif (numbers =6.3) then
plot3d([2*cos(t), sin(t), -2+2*cos(t)+sin(t)], t=0..2*Pi, s= 0..1, grid=[30,2],
axes= NORMAL) ;

elif (numbers = 6.4) then
plot3d([cos(t), cos(t) -1 +2*sin(t), sin(t)], t=0..2*Pi, s= 0..1, grid=[30,2],
axes= NORMAL) ;

elif (numbers = 6.5) then
picture1 := plot3d([t, t^3, t], t = -2..2, s=0..1, grid=[30,2]) :
picture2 := plot3d( [ 1 +t, 1+3*t, 1+t], t=0..1, s=0..1, grid=[2,2]) :
plots[display3d]( {picture1, picture2} ) ;

elif (numbers = 6.6) then
plot3d([t, abs(t), sin(t)], t=-2*Pi..2*Pi, s= 0..1, grid=[30,2], axes= NORMAL)
;

elif (numbers = 6.8) then
curve1 := plot3d( [1 +5*cos(t)*3/5+5*4/13*sin(t), 1 + 5*cos(t)*4/5
-3/13*5*sin(t), -1 + 12/13*5*sin(t)], t=0..2*Pi, s=0..1, grid =[60,2]) :
graph1 := plot3d( ( 48*x - 36*y -37 )/ 25, x=-10..10, y=-10..10, grid=[4,4]) :
with(plots) :  display3d( {curve1, graph1} ) ;  

elif (numbers = 7.1) then
plot3d( (5 - 2*x+3*y )/2, x=-2..2, y=-2..2, axes= NORMAL) ;

elif (numbers =7.2) then
plot3d({sqrt((1-x^2/4-y^2)*9), -sqrt((1-x^2/4-y^2)*9)}, x=-2..2,y=-2..2) ;

elif (numbers = 7.3) then
plot3d({x+sqrt(1-y^2), x-sqrt(1-y^2)}, x=-2..2, y=-2..2) ;  

elif (numbers = 7.4) then
plot3d( { sqrt(x^2-1+y^2/4), -sqrt(x^2-1+y^2/4), sqrt(x^2+y^2/4+1),
-sqrt(x^2+1+y^2/4)}, x=-2..2, y=-2..2) ; 

elif (numbers =7.5) then
plot3d( x^2*y^2*(x^2-y^2), x=-2..2, y=-2..2, orientation=[52,70]) ;

elif( numbers = 10.1) then
plot3d (x^3 +y ,x=-2..2, y=-2..2) ; 

elif(numbers = 10.2) then
plot3d( [2*cos(t), v, 3*sin(t)] , t=0..2*Pi, v=-2..2) ;

elif(numbers = 10.3) then
plot3d( [x, cos(x*z) +1, z], x=-2..2, z=-2..2) ; 

elif(numbers = 10.4) then
plot3d( [2*sin(v)*cos(u),  3*sin(v)*sin(u), 5*cos(v)], v=0..Pi, u=0..2*Pi) ; 

elif(numbers =10.5) then
plot3d( { [sinh(u)*cos(v), sinh(u)*sin(v), cosh(u)], [sinh(u)*cos(v),
sinh(u)*sin(v),  -cosh(u)]}, v=0..2*Pi, u=-2..2) ;

elif (numbers = 10.6)  then
plot3d( [ (4+sin(v))*sin(2*u), (4+sin(v))*cos(2*u), 2*u+cos(v)], u= 0..2*Pi,
v=0..2*Pi, grid = [35,20]) ;  

elif(numbers = 10.7) then
plot3d( [ u*cos(v), u*sin(v), 2*u*cos(v)+3*u*sin(v)], u=0..1, v=0..2*Pi, grid =
[5, 40], orientation = [30, 70]) ;

elif(numbers = 10.8) then
plot3d( [u*cos(v), u*sin(v), sqrt( (4 - (u*cos(v))^2 - (u*sin(v))^2)/2)],
u=0..1, v=0..2*Pi, grid = [5, 40] );

elif(numbers = 11.1) then
plot3d( {[r*cos(t), r*sin(t), r], [r*cos(t), r*sin(t), -r] }, r=0..2,
t=0..2*Pi, grid=[10,40]);

elif(numbers = 11.2) then
plot3d( { [ r*cos(t), r*sin(t), sqrt(4-4*r^2)], [r*cos(t), r*sin(t),
-sqrt(4-4*r^2)] }, r=0..1, t=0..2*Pi, grid=[30,30] ) ;

elif(numbers = 11.3) then
plot3d( { [r*cos(t), r*sin(t), r^2], [r*cos(t), r*sin(t), 4 - r^2]},
r=0..sqrt(2), t=0..2*Pi, grid=[30, 20]) ;

elif(numbers=11.4) then
plot3d( [ r*cos(Pi/2), r*sin(Pi/2), z], r=0..8, z=0..2, axes =NORMAL) ;

elif(numbers=11.5) then
plot3d([r*cos(t), r*sin(t), 1], r=0..1, t=0..2*Pi, grid=[10,30], axes=NORMAL);
elif(numbers=11.6) then
plot3d([1*cos(t), 1*sin(t), z],  t=Pi/2..2*Pi, z=0..1, axes= NORMAL ) ;

elif(numbers=11.7) then
plot3d([r*cos(t), r*sin(t), 0], r=1..2, t=0..Pi, axes=NORMAL, scaling =
CONSTRAINED) ;

elif(numbers=11.8)  then
plot3d( [r*cos(t), r*sin(t), cos(r)^2], r=0..2*Pi, t=0..2*Pi ) ;

elif(numbers=11.9) then
plot3d( [cos(t)^2*cos(t), cos(t)^2*sin(t), z], t =0..2*Pi, z=0..1) ; 

elif(numbers=11.11) then
plot3d( [2*cos(t), 2*sin(t), z], t=0..3*Pi/2, z=-1..1, axes = NORMAL, scaling =
CONSTRAINED) ; 

elif(numbers=11.12) then
plot3d( [r*cos(Pi/2), r*sin(Pi/2), z], r=1..2, z=0..2, axes=NORMAL, scaling =
CONSTRAINED) ;

elif(numbers=11.3) then
plot3d([ r*cos(t), r*sin(t), 2], r=0..1, t=Pi/2..Pi, axes=NORMAL, scaling=
CONSTRAINED) ;

elif(numbers=12.1) then
plot3d([cos(t)*sin(p)*cos(t), cos(t)*sin(p)*sin(t), cos(t)*cos(p)], t=0..2*Pi,
p=0..Pi) ; 

elif(numbers=12.2) then
plot3d( [ (5*cos(p)^2 -1)/2 *sin(p)*cos(t), (5*cos(p)^2 -1)/2 *sin(p)*sin(t),
(5*cos(p)^2 -1)/2*cos(p)], t=0..2*Pi, p=0..Pi) ; 

elif(numbers=12.3) then
plot3d( [ cos(t)^2, cos(t)*sin(t), cos(t)*cos(p)/sin(p)], t=0..2*Pi, p=0..Pi,
grid=[30,10]) ;

elif(numbers=12.4) then
plot3d([ sqrt(1/(sin(p)^2 - cos(p)^2))*sin(p)*cos(t), sqrt(1/(sin(p)^2 -
cos(p)^2))*sin(p)*sin(t), sqrt(1/(sin(p)^2 - cos(p)^2))*cos(p)], t=0..2*Pi,
p=0..Pi) ;

elif(numbers=12.5) then
plot3d([
sqrt(1/(sin(p)^2*cos(t)^2+2*sin(p)^2*sin(t)^2+4*cos(p)^2))*sin(p)*cos(t),
sqrt(1/(sin(p)^2*cos(t)^2+2*sin(p)^2*sin(t)^2+4*cos(p)^2))*sin(p)*sin(t),
sqrt(1/(sin(p)^2*cos(t)^2+2*sin(p)^2*sin(t)^2+4*cos(p)^2))*cos(p)], t=0..2*Pi,
p=0..Pi) ;

elif(numbers=12.6) then
plot3d( [ abs(1/(sin(p)*cos(t)+2*sin(p)*sin(t)+cos(p)))*sin(p)*cos(t),
abs(1/(sin(p)*cos(t)+2*sin(p)*sin(t)+cos(p)))*sin(p)*sin(t),
abs(1/(sin(p)*cos(t)+2*sin(p)*sin(t)+cos(p)))*cos(p)], t=0..2*Pi, p=0..Pi,
grid=[30,30]) ;

elif(numbers=12.7) then
plot3d( [2*sin(p)*cos(t), 2*sin(p)*sin(t), 2*cos(p)], p=0..Pi/2, t=0..3*Pi/2,
axes =NORMAL) ; 

elif(numbers=12.8) then
plot3d( [r*sin(Pi/6)*cos(t), r*sin(Pi/6)*sin(t), r*cos(Pi/6)], r=0..1, t=0..Pi,
axes=NORMAL, scaling= CONSTRAINED) ; 

elif(numbers=12.9) then
plot3d( [r*sin(p)* cos(Pi/2), r*sin(p)*sin(Pi/2), r*cos(p)], r=0..2,
p=Pi/4..3*Pi/4, axes=NORMAL) ; 

elif(numbers=12.11) then
plot3d( [2*sin(p)*cos(t), 2*sin(p)*sin(t), 2*cos(p)], t=0..Pi, p=Pi/2..Pi, axes
= NORMAL, scaling = CONSTRAINED) ;

elif(numbers=12.12) then
plot3d( [ r*sin(3*Pi/4)*cos(t), r*sin(3*Pi/4)*sin(t), r*cos(3*Pi/4)], r=0..3,
t=Pi/2..3*Pi/2, axes=NORMAL, scaling = CONSTRAINED) ;

elif(numbers=12.13) then
plot3d( [ r*sin(p)*cos(Pi), r*sin(p)*sin(Pi), r*cos(p)], r=1..3, p=0..3*Pi/4,
axes=NORMAL, scaling = CONSTRAINED) ;  

elif (numbers = 13.1)  then
plot3d([x, y, x^2+4*y^2-sin(x*y) - 2*y],  x=-1..2, y=-1..2, axes = NORMAL) ; 



elif (numbers = 13.4)  then
plot( (x+0.5)*(x-1/2)*(x-1), x=-1..1.5) ; 

elif (numbers =14.2)  then
f := x^2 +x*y :
v := [diff(f,x), diff(f,y)] : with(mvcal) :
plotvectorfield(v, x=-1..1, y=-1..1, 1) ; 

elif(numbers = 15.1) then
picture := plot3d( 5*sqrt( 2 - x^2/4- y^2/9), x=-2..2, y=2..4) :
tangent := plot3d( -5/3*y+10, x=-1..1, y=2..4, grid = [5, 5]) :
plots[display3d]( {picture, tangent} ) ;  

elif(numbers = 15.2) then
picture := plot3d( x*y-x^2*y+1, x=0..2, y=0..2, scaling= CONSTRAINED) :
normal := plot3d( [1-1*t, 1+0*t, 1-t], t=-5..5, s=0..1, grid=[2,2]) :
plots[display3d]( {picture, normal} ) ;

elif(numbers = 15.3) then
picture := plot3d([(4+cos(v))*cos(u), (4+cos(v))*sin(u), sin(v)], u=0..2*Pi,
v=0..2*Pi)  :  
tangent := plot3d( 1, x=-1..1, y=3..5, grid=[5,5]) :
plots[display3d]( {picture,tangent}) ;

elif(numbers = 15.4) then
plot3d( {x^2+2*x*y+2*y^2, 5}, x=-1.5..1.5, y=-1.5..1.5) ;

elif(numbers=15.5) then
plot3d( {[t, -1/2*t+1/2*(-t^2+10)^(0.5), 5], [1+4*s, 1+6*s,5-s],
[1-6*s,1+4*s,5]}, t=-2..2, s=0..1, grid=[30,2], scaling= CONSTRAINED) ;  

elif(numbers=15.6) then
picture := plot3d( [2*u^2, u*cos(v), u*sin(v)], u=0..1.5, v=-Pi/2..Pi/2, grid=
[5,30]):   line := plot3d( [-1/2+t, 13-4*t, 13/2-2*t], t=0..3, s=0..1,
grid=[2,2]) :  plots[display3d]({picture, line} ) ;

elif(numbers = 16.4) then
f :=y^4 - 2*y^2 - x^2 +x;  
plotvectorfield( [ diff(f, x)/4, diff(f,y)/4], x=0..1, y=-1.2..1.2, 0.2) ;

elif(numbers=17.3) then
f := x^2-2*y ;  
picture := plotvectorfield( [diff(f,x)/4, diff(f,y)/4], x=-2..2, y=-2..2, 0.5)
:
circle := plot( [2*cos(t), 2*sin(t), t=0..2*Pi] ) :  plots[display]( { picture,
circle} ) ;

elif(numbers=17.5) then
plot( [1/3*(2-cos(2*(t+2*Pi/3)))*cos(t),
1/4*(2-cos(2*t))*sin(t), t=0..2*Pi ], x=-1..1, y=-1..1) ; 

elif(numbers=18.4) then
dzdydxplot(z=0..4+x^2+y^2, y=1..x^2, x=2..3) ; 

elif(numbers = 18.5) then
dzdxdyplot(z=0..x^2+2*y^2, x=-sqrt(1-(y-1)^2)+1..sqrt(1-(y-1)^2)+1, y=0.5..1.5)
;

elif (numbers = 20.1) then
plot([(1-cos(theta))*cos(theta), (1-cos(theta))*sin(theta), theta=0..2*Pi]);  

elif (numbers = 20.2)  then
plot({[(1-cos(theta))*cos(theta), (1-cos(theta))*sin(theta), theta = 0..2*Pi] ,
 [(1+cos(theta))*cos(theta),  (1+cos(theta))*sin(theta), theta =0..2*Pi]});  

elif (numbers = 20.3)  then
plot ({[sqrt(2*cos(2*theta))*cos(theta),
sqrt(2*cos(2*theta))*sin(theta),theta=-Pi..3*Pi], [cos(theta), sin(theta),
theta=0..2*Pi] }, x=-2..2,y=-2..2);  

elif(numbers=20.5) then
drdtplot( r=theta^2..(theta+2*Pi)^2, theta=0..2*Pi ) ;  

elif(numbers=21.1) then
dzdydxplot( z=0..(1-x-2*y)/3, y=0..(1-x)/2, x=0..1) ;

elif(numbers=21.2) then
dxdzdyplot( x=0..2-z-y,  z=sqrt(y)..y^2, y=0..1) ;

elif(numbers=21.3) then
dydzdxplot( y=0..2-2*z, z=0..(2-x)/2, x=0..2) ;

elif(numbers=21.4) then
dydzdxplot( y=0..sqrt(4-x^2-z), z=0..4-x^2, x=0..2) ; 

elif(numbers=21.5) then
dxdydzplot( x=0..(3-z)/3, y=z/3..(z-6)/(-3), z=0..3) ; 

elif(numbers=21.6) then
dzdxdyplot( z= -sqrt( (1.5)^2-x^2-(y-1)^2) +6..4-x^2-y^2, x=-1..1, y=1..2) ;

elif(numbers=21.7) then
dydzdxplot( y=0..sqrt(1-z^2), z=0..sqrt(1-x^2), x=0..1) ; 

elif(numbers=21.8) then
dydxdzplot( y=0..3, x=z^2..2-z^2, z=-1..1); 

elif(numbers=21.9) then
dydzdxplot( y=(3*z-x)/2..2, z=x^2..(x+4)/3, x=-1..4/3) ; 

elif(numbers=21.12) then
dydxdzplot( y=(z-4)/2..(-z+4)/2, x=(z-4)/2..(-z+4)/2, z=0..2) ; 

elif(numbers=21.13) then
dydzdxplot( y= x^2+z^2..2-x^2-z^2, z=-sqrt(1-x^2)..sqrt(1-x^2), x=-1..1 ) ; 

elif(numbers = 21.14) then
sphere := plot3d( [ 2*sin(p)*cos(t), 2*sin(p)*sin(t), 2*cos(p)], p=0..Pi,
t=0..2*Pi)  :
plane := plot3d( 1/2*y+2, x=-2..2, y=-2..2) : 
plots[display3d]( {sphere, plane}, scaling= CONSTRAINED ) ;

elif(numbers = 21.15) then
dzdxdyplot( z=1/2*y+2..sqrt(4-x^2-y^2),
x=-1/2*(-5*y^2-8*y)^(1/2)..1/2*(-5*y^2-8*y)^(1/2), y=-8/5..0) ;  

elif (numbers = 21.13)  then
graph1 := plot3d( [sqrt(2)*sin(p)*cos(t), sqrt(2)*sin(p)*sin(t),
sqrt(2)*cos(p)], p=0..Pi, t=0..2*Pi) :
graph2 := plot3d( y^2, x=-2..2, y=-2..2) : with(plots):
display3d({graph1,graph2} ) ;  

elif(numbers=22.1) then
stair := plot3d([r*cos(t), r*sin(t), t], t=0..Pi/2, r=0..1) :
picture := plot3d({[cos(t*Pi/2), sin(t*Pi/2), s*t*Pi/2], [0,t,s*Pi/2]}, t=0..1,
s=0..1) :   plots[display3d] ({stair, picture} ) ;

elif(numbers=22.2) then
graph1 := plot3d({[4*s*cos(t)*cos(t), 4*s*cos(t)*sin(t), -sqrt(16 -
(4*s*cos(t))^2)], [4*s*cos(t)*cos(t), 4*s*cos(t)*sin(t), sqrt(16 -
(4*s*cos(t))^2)]}, s=0..1, t=0..Pi/2, grid=[20,20]) :  
graph2 := plot3d( [4*cos(t)*cos(t), 4*cos(t)*sin(t), -s*sqrt(16 -
(4*cos(t))^2)+(1-s)* sqrt(16 - (4*cos(t))^2)], t=0..Pi/2, s=0..1, grid=[20,
10]) :  graph3 := plot3d({[4*t, 0, -s*sqrt(16 - (4*t)^2)+(1-s)* sqrt(16 -
(4*t)^2)], [0, 4*t, -s*sqrt(16 - (4*t)^2)+(1-s)* sqrt(16 - (4*t)^2)]}, t=0..1,
s=0..1, grid=[5,10]) :
plots[display3d]({graph1, graph2}, scaling= CONSTRAINED) ;  

elif(numbers=22.3) then
graph1 := plot3d([r*cos(t), r*sin(t),sqrt(4-3*r^2)], r=0..1, t =0..Pi/2) :
graph2 := plot3d( [r*cos(t), r*sin(t), r^2], r=0..1, t=0..Pi/2) : 
graph3 := plot3d( {[s, 0, t*sqrt(4-3*s^2)+(1-t)*s^2], [0, s,
t*sqrt(4-3*s^2)+(1-t)*s^2]}, s=0..1, t=0..1, grid=[5,10]) :
plots[display3d] ({graph1, graph2, graph3} ,scaling= CONSTRAINED) ;

elif(numbers=22.4) then
graph1 := plot3d( {[ 1*cos(t), 1*sin(t), z] , [2*cos(t)*cos(t),
2*cos(t)*sin(t),z] }, t=-Pi/3..Pi/3, z=0..2):
graph2 := plot3d( {[ (s*2*cos(t)+(1-s))*cos(t), (s*2*cos(t)+(1-s))*sin(t), 2],
[ (s*2*cos(t)+(1-s))*cos(t), (s*2*cos(t)+(1-s))*sin(t), 0]}, t=-Pi/3..Pi/3,
s=0..1, grid=[20,5]) :
plots[display3d]({graph1, graph2} ) ;

elif(numbers=22.5) then
graph1 := plot3d( {[ 1*cos(t), 1*sin(t), z] , [2*sin(t)*cos(t),
2*sin(t)*sin(t),z] ,[-2*sin(t)*cos(t), -2*sin(t)*sin(t),z]  }, t=-Pi/6..Pi/6,
z=0..2):  graph2 := plot3d([cos(t), sin(t) ,z], t=Pi-Pi/6..Pi+Pi/6, z=0..2) : 
plots[display3d]({graph1, graph2} ) ;

elif(numbers=22.7) then
plot3d( [r*cos(t), r*sin(t), 3 - (r^2-0.5^2)*(r^2-0.7^2)*(r^2-1)], r=0..1,
t=0..2*Pi) ;

elif(numbers=23.1) then
sphere := plot3d([ 4*cos(p)*sin(p)*cos(t), 4*cos(p)*sin(p)*sin(t),
4*cos(p)*cos(p)], t=0..2*Pi, p=0..Pi/6) :  cone := plot3d([r*sin(Pi/6)*cos(t),
r*sin(Pi/6)*sin(t), r*cos(Pi/6)], r=0..4*cos(Pi/6), t=0..2*Pi) : 
plots[display3d]( {sphere, cone}) ; 

elif(numbers=23.2) then
cone := plot3d([r*sin(Pi/4)*cos(t), r*sin(Pi/4)*sin(t), r*cos(Pi/4)],
r=0..2/cos(Pi/4), t=0..2*Pi ) :  cover := plot3d([r*cos(t), r*sin(t),
2],r=0..2, t=0..2*Pi)  :  plots[display3d]({cone, cover}) ;

elif(numbers=23.3) then
sphere := plot3d([4*cos(p)*sin(p)*cos(t), 4*cos(p)*sin(p)*sin(t),
4*cos(p)*cos(p)], p=arccos(3/4)..Pi/2, t=0..2*Pi):  sphere2:=
plot3d([3*sin(p)*cos(t), 3*sin(p)*sin(t), 3*cos(p)], p=0..arccos(3/4),
t=0..2*Pi) :  plots[display3d]({sphere, sphere2}, scaling=CONSTRAINED) ;

elif(numbers=23.4) then
sphere := plot3d([4*cos(p)*sin(p)*cos(t), 4*cos(p)*sin(p)*sin(t),
4*cos(p)*cos(p)], p=0..arccos(1/2), t=0..2*Pi):  cover := plot3d([r*cos(t),
r*sin(t), 1], r=0..sqrt(3), t=0..2*Pi) : plots[display3d]({sphere, cover},
scaling=CONSTRAINED) ;

elif(numbers=23.5) then
plot3d( [(1-cos(p))*sin(p)*cos(t), (1-cos(p))*sin(p)*sin(t),
(1-cos(p))*cos(p)], p=0..Pi, t=0..2*Pi, scaling=CONSTRAINED) ;

elif(numbers = 23.7) then
plot3d([ 0.7*cos(3*p)*sin(p)*cos(t),  0.7*cos(3*p)*sin(p)*sin(t), 
0.7*cos(3*p)*cos(p)], p=0..Pi/6, t=0..2*Pi, scaling= CONSTRAINED) ;

elif(numbers=24.6) then
plotvectorfield( [cos(x*y), sin(y)+1/2], x=-0.5..2, y=-0.5..2, 0.5);

elif(numbers=25.1) then
plotvectorfield([y, -x], x=-2..2, y=-2..2, 0.5) ;

elif(numbers=25.2) then
plotvectorfield([0.3, 0], x=-2..2, y=-2..2, 0.5) ;

elif(numbers=25.3) then
plotvectorfield([0, (x+y)/4], x=-2..2, y=-2..2, 0.5) ;

elif(numbers=28.3) then
plot3d( [u^2*cos(t), u, u^2*sin(t)], t=0..2*Pi, u=0..2) ;

elif(numbers=28.4) then
plot3d( [ (4+cos(v))*cos(u), (4+cos(v))*sin(u), 2*u+sin(v)], u=0..2*Pi,
v=0..2*Pi) ; 

elif (numbers = 4.9)  then
plot3d( { [103*t+(1-t)*733, 51*t+(-369)*(1-t), 1254*t], [103*t+(1-t)*418,
51*t+(471)*(1-t), 1254*t], [103,51,1254*t],[0,0,0],[800*s,850*t-400,0]},
t=0..1, s=0..1, grid=[5,5],axes=NORMAL, orientation=[41,63]) ; else

fi; 
end :

mvcal[eg3] := proc(numbers);
if (numbers =4.1) then
vline := [0,0,1,2]; wline := [0,0,3,3]; segment := [1,2,3/2,3/2];vhead1:= [1,2,
1, 1.8]; vhead2 := [1, 2, 0.8,1.9]; whead1 := [1, 1, 0.9, 0.8]; whead2 :=
[1,1,0.8,0.9];
 plot({vline,wline,segment, vhead1, vhead2, whead1, whead2});  fi; end :

mvcal[eg5] := proc(numbers) ;
if (numbers =2.2) then
with(plots) :
vector1 := plot3d([0-s, 1, Pi/2+s], t=0..1, s=0..1, grid=[2,2]) : 
vector2 := plot3d([1/sqrt(2) -1/sqrt(2)*s, 1/sqrt(2) + 1/sqrt(2)*s, Pi/4 + s],
t=0..1, s=0..1, grid=[2,2]) :
helix := plot3d([ cos(t), sin(t), t], t=0..2*Pi, s=0..1, grid =[40,2]) :
display3d({ vector1, vector2, helix} ) ;   fi; end :

mvcal[eg6] := proc(numbers) ;
if (numbers = 1.7) then
plane := plot3d( 1 -2*x +y, x=-2..2, y=-2..2) :
cylinder := plot3d( [2*cos(t), 2*sin(t), s], t=0..2*Pi, s=-7..7) :
with(plots) :
display3d( {plane, cylinder} );  fi; end :

mvcal[eg7] := proc(numbers) ;  
if (numbers = 3.1) then
plot3d( [8*cos(t)*sin(p), 8*sin(t)*sin(p), 8*cos(p)], t=0..2*Pi, p=0..Pi,
scaling = CONSTRAINED, axes = NORMAL);

elif (numbers = 3.2) then
plot3d( [ 1*sin(p)*cos(t), 2*sin(p)*sin(t), 3*cos(p)], p=0..Pi, t=0..2*Pi,
grid= [30,30] , scaling = CONSTRAINED,axes= NORMAL ) ;  

elif (numbers = 3.3)  then
plot3d(  {[x, y, -x^2-y^2], [0, 0,0]}, x=-2..2, y=-2..2, axes = NORMAL) ;  

elif (numbers = 3.4) then
plot3d(  {[x, y, x^2+y^2 + 3], [0, 0,0]}, x=-2..2, y=-2..2, axes = NORMAL) ;  

elif (numbers = 3.5)  then
plot3d( [ r*sin(Pi/4)*cos(t), 2*r*sin(Pi/4)*sin(t), 3*r*cos(Pi/4)], r=-4..4,
t=0..2*Pi, grid= [30,30] , scaling = CONSTRAINED,axes= NORMAL ) ;  else 

fi;  end:

mvcal[eg8] := proc (numbers) ;
if (numbers = 1.1) then
plot3d({ [0,0,0], [sqrt(2)*cos(t), sqrt(2)*sin(t), 2], [2*cos(t), 2*sin(t),
4]}, t=0..2*Pi, s=0..1, grid= [40,2], title = `contour curves`) ;  else

if (numbers = 1.2)  then
plot( {[0,0, t=0..2*Pi], [sqrt(2)*cos(t), sqrt(2)*sin(t), t=0..2*Pi],
[2*cos(t), 2*sin(t), t=0..2*Pi]}, title = `level curves`) ;  else

if (numbers = 1.5)  then
plot3d( { [r*cosh(t), r*sinh(t), (cosh(t))^2 + (sinh(t))^2], [r*sinh(t),
r*cosh(t), -(cosh(t))^2 - (sinh(t))^2 ]}, r=-1..1, t=-1..1, axes = NORMAL ) ; 

fi; fi; fi; end :

mvcal[eg9] := proc(numbers) ;
if (numbers = 2.1) then
plot3d([cosh(v)*cos(u),cosh(v)*sin(u),sinh(v)],u=0..2*Pi,v=-1..1); else

if (numbers = 2.2) then
plot3d( [sinh(v)*cos(u), sinh(v)*sin(u), sinh(v)] , u=0..2*Pi, v=-1..1);  else

if (numbers = 2.3) then
plot3d( {[sinh(v)*cos(u), sinh(v)*sin(u), cosh(v)], [sinh(v)*cos(u),
sinh(v)*sin(u), -cosh(v)] }, u=0..2*Pi, v=-1..1) ;  else

if (numbers = 2.4) then
plot3d( {[sinh(v)*cos(u), sinh(v)*sin(u), cosh(v)], [sinh(v)*cos(u),
sinh(v)*sin(u), -cosh(v)], [cosh(v)*cos(u),cosh(v)*sin(u), sinh(v)],
[sinh(v)*cos(u), sinh(v)*sin(u), sinh(v)] }, u=0..2*Pi, v=-1..1);  
fi;fi; fi; fi; end :

mvcal[eg11] := proc (numbers) ;
if (numbers = 1.8) then
fig6 := plot3d([r*cos(t), r*sin(t), r^2], r=0..2, t=0..2*Pi, grid=[5,40]);
fig7 := plot3d([t, sqrt(2*Pi*t- t^2), z], t=0..2*Pi, z=-2..2, grid=[40,6]);
fig8 := plot3d([sqrt(2*Pi*t-t^2)*cos(t),  sqrt(2*Pi*t-t^2)*sin(t), z],
t=0..2*Pi, z=-2..2, grid=[60,5]); with(plots) :  fi;
end :

mvcal[eg12] := proc (numbers) ;
if (numbers=7) then
sphere := plot3d( [2*sin(phi)*cos(theta), 2*sin(phi)*sin(theta),2*cos(phi)],
phi=0..Pi, theta=0..2*Pi) :
cone := plot3d([rho*sin(Pi/4)*cos(theta), rho*sin(Pi/4)*sin(theta),
rho*cos(Pi/4)], rho =0..4, theta =0..2*Pi) :
with(plots) :
display3d( {sphere, cone});  fi ; end:

mvcal[eg13] := proc(numbers);
if (numbers = 1.4) then 
graph1 := plot3d ( [ x, 2, 4/(2+cos(x))], x=Pi -1.. Pi +1, y=1..3) :
 tangent1 := plot3d ([Pi +t, 2, 4+t*0], t=0..1, x=Pi-1..Pi+1, grid=[2,2]) :
graph2 := plot3d( [Pi, y, 2*y/(y+cos(Pi) )], y=1..3, x=Pi-1..Pi+1) :
 tangent2 := plot3d([Pi, 2+t, 4+t*(-2)], t=0..1, y=1..3, grid= [2, 2]) :
with(plots) :  
display3d( {graph1, graph2, tangent1, tangent2}, axes=NORMAL) ;  fi; end :

mvcal[eg14] := proc(numbers) ;
if (numbers = 2.1) then
plot({[1,0,1,1],[0,1,-1,1],[-1,0,-1,-1],[0,-1,1,-1],
[1.5,1.5,0,3],[-1.5,1.5,-3,0],[-1.5,-1.5,0,-3],[1.5,-1.5,3,0],
[1,1,1.2,0.8],[1,1,.8,.8],[-1,1,-.8,.8],
[-1,1,-.8,1.2],[-1,-1,-1.2,-.8],[-1,-1,-.8,-.8],[1,-1,.8,-1.2],[1,-1,.8,-.8],
[0,3,.2,2.6],[0,3,.2,3],[-3,0,-2.8,.4],[-3,0,-2.8,0.1],[3,0,2.8,-0.4],
[3,0,2.8,0.1],[0,-3,-0.2,-3],[0,-3,-0.2,-2.6]}) ;  else

 if (numbers = 3.1) then 
graph := plot3d(x*y+y^3, x=0.5..1.5, y=0.5..1.5, orientation=[40,80]) :
with(plots) ;
 graph3 :=  plot3d([1 + t/sqrt(5), 1+2*t/sqrt(5), s*((1 + 
t/sqrt(5))*(1+2*t/sqrt(5)) + (1+2*t/sqrt(5))^3)], s=0..1, t=0..0.5) :
 display3d( { graph, graph3} );   else

if (numbers = 3.2) then
graph1 := plot3d(sin(x*y)+1, x=-2..0, y=2..4, orientation=[74,79]) :
  direction1 := plot3d([-1+2/sqrt(13)*t, Pi + 3/sqrt(13)*t,
s*(sin((-1+2/sqrt(13)*t)*(Pi + 3/sqrt(13)*t))+1)], s=0..1, t=0..0.5) :
  direction2 := plot3d( [-1 + 1/sqrt(50)*t, Pi +7/sqrt(50)*t, s* (sin( (-1 +
1/sqrt(50)*t)*(Pi +7/sqrt(50)*t))+1)], s=0..1, t=0..0.5) :
  display3d( {graph1,direction1, direction2}) ;   else

if (numbers = 5) then
levelcurves :=  levelcurve( f, [0,1,2,3,4], x=-2..2, y=-2..2) :
gradient :=  plotvectorfield( [ diff(f, x) , diff(f, y)], x=-2..2, y=-2..2,
0.5) ;
with (plots) :  display({levelcurves}) ;
display({gradient}) ;  else

if (numbers = 5.1) then
curves :=  levelcurve( f, [0,1,4,9], x=-3..3, y=-3..3) :
gradient :=  plotvectorfield( [ diff(f, x) , diff(f, y)], x=-2..2, y=-2..2,
0.5):
with (plots) :   display({curves}) ;  else

if (numbers = 5.2) then
curves := levelcurve(f, [0,1,2,3, -1, -3, -7,-12], x=-3..3, y=-3..3) :
gradient := plotvectorfield( [diff(f,x)/10, diff(f,y)/10], x=-2..2, y=-2..2, 1)
:
with (plots) :   display({curves}) ;  else

if (numbers =5.3) then
curves := levelcurve( f, [0,-5, -10,-20, 10, 20,40], x=-3..3, y=-3..3) :
gradient := plotvectorfield([diff(f, x)/20 , diff(f,y)/20], x=-2..2, y=-1..1,
0.5):   plots[display]({curves, gradient}) ; 

fi ; fi; fi; fi; fi; fi; fi; end : 

mvcal[eg19] := proc ( numbers) ;
if (numbers = 1.1)  then
plot3d( { [v*cos(u), v*sin(u), 6 -2*(v*cos(u))^2-(v*sin(u))^2 ], [v*cos(u),
v*sin(u),
 2*(v*sin(u))^2 + (v*cos(u))^2], [cos(u), sin(u), 
v*(6 -2*(1*cos(u))^2-(1*sin(u))^2)+ (1-v)*(2*(1*sin(u))^2 + (1*cos(u))^2)]},
v=0..1, u=0..2*Pi, grid=[10,30]);  else

if (numbers = 1.2) then
plot3d( { [v*cos(u), v*sin(u), 6 -2*(v*cos(u))^2-(v*sin(u))^2 ], [cos(u),
sin(u), 
v*(6 -2*(1*cos(u))^2-(1*sin(u))^2)] }, v=0..1, u=0..2*Pi, grid=[10,30]);  else

if (numbers = 1.3) then 
plot3d( { [v*cos(u), v*sin(u), 2*(v*sin(u))^2 + (v*cos(u))^2], [cos(u), sin(u),

(1-v)*(2*(1*sin(u))^2 + (1*cos(u))^2)] }, v=0..1, u=0..2*Pi, grid=[10,30]); 
else

if (numbers =1.4) then
plot3d({[u,v,u+v],[u,v,u+7*v],[0,u,u+6*u*v],[1,u,1+u+6*u*v],
  [u,1,u+1+6*v]},u=0..1,v=0..1, grid=[10,10],
  orientation=[15,60]); 

fi; fi ; fi; fi; end :

mvcal[eg22] := proc (numbers); 
if (numbers =1.2) then 
graph1 := plot3d ( { [2*cos(theta), 2*sin(theta), 2*r], [ r*cos(theta), 
r*sin(theta), r^2],[r*cos(theta), r*sin(theta), 0]}, theta = Pi/2..2*Pi, r =
0..2 ):
graph2 := plot3d({ [u,0,s*u^2], [0,u,s*u^2]}, u=0..2, s=0..1) :
 with(plots);  display3d ( { graph1, graph2});  else

if (numbers = 1.4) then
graph1 := plot3d([cosh(v)*cos(theta),cosh(v)*sin(theta),sinh(v)],  
theta=Pi/2..2*Pi,v=ln(sqrt(2)-1)..ln(sqrt(2)+1)):
graph2 := plot3d([sqrt(2)*cos(theta),
 sqrt(2)*sin(theta), v], theta=Pi/2..2*Pi,v=-1..1) :
 graph3 := plot3d({[0,u, s*(-sqrt(u^2-1)) +(1-s)*sqrt(u^2-1)], [u,0,s*(-
sqrt(u^2-1)) +(1-s)*sqrt(u^2-1)] } , u=1..sqrt(2), s=0..1) :
with(plots):
display3d( { graph1, graph2, graph3} );  else

if (numbers =1.6) then
graph1 := plot3d( [(cos(2*theta))^2*cos(theta), (cos(2*theta))^2*sin(theta),
z], theta =0..2*Pi, z=0..2 , grid = [60, 5]):
graph2 := plot3d( {[r*(cos(2*theta))^2*cos(theta),
r*(cos(2*theta))^2*sin(theta),2], [r*(cos(2*theta))^2*cos(theta),
r*(cos(2*theta))^2*sin(theta), 0] }, theta=0..2*Pi, r=0..1, grid =[50,10]) :
with(plots);  display3d({ graph1, graph2});  

fi ; fi; fi; end :

mvcal[eg23] := proc( numbers) 
if (numbers = 1.1) then
graph1 := plot3d([2*sin(v)*cos(u),2*sin(v)*sin(u),2*cos(v)],u 
=Pi/2..2*Pi, v = Pi/6..5*Pi/6):
graph2 := plot3d([cos(u), sin(u), z], u=Pi/2..2*Pi, z=-sqrt(3)..sqrt(3)):
graph3 := plot3d({[0,v,t*(-sqrt(4-v^2))+(1-t)*sqrt(4-v^2)], [v,0,t*(-
sqrt(4-v^2))+(1-t)*sqrt(4-v^2)]}, v=1..2, t=0..1) :
 with(plots):
display3d({graph1, graph2, graph3}); else

if (numbers = 1.4) then
graph1 := plot3d(  [2*cos(phi)*sin(phi)*cos(theta), 
2*cos(phi)*sin(phi)*sin(theta),2*cos(phi)*cos(phi) ], phi=Pi/3..Pi/2, theta 
=Pi/2..2*Pi ) :
graph2 := plot3d([ sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)]  , 
phi=0..Pi/3, theta =Pi/2..2*Pi ) :
graph3 := plot3d( { [0,u, t*sqrt(1-u^2)+(1-t)*(1- sqrt(1-u^2))], [u,0,t*sqrt(1-
u^2) + (1-t)*(1-sqrt(1-u^2))]}, u=0..sqrt(3/4), t=0..1) :
with(plots) :
display3d({graph1, graph2, graph3});  
fi; fi; end:

 
mvcal[eg27] := proc( numbers) ;
if (numbers = 1.2) then
plot3d( {[t, s*(1-2*t), 2/(1+t)], [t, s*(1-2*t),0]}, t=0..1/2, s=0..1) ;  else

if (numbers = 2.1) then
plotsvf3d([u, v, u+v-4], [x, 2*y, z], u=0..1, v=0..2);  else

if (numbers = 2.2) then
plotsvf3d([1/sqrt(2)*cos(u), v, sin(u)], [y*z,x+y,z^2], u=0..Pi, v=-1..1); else

if (numbers =2.3) then
plotsvf3d([u,v,0], [x^2,y+z,y], u=-1..1,v=-1..1);  

fi;  fi; fi; fi; end :

mvcal[eg28] := proc( numbers);
if (numbers = 3.1) then plot3d([s*(2-cos(4*t))*cos(t), s*(2-cos(4*t))*sin(t),
s^2], t=0..2*Pi, 
s=0..2, grid= [50,10]);  else

if (numbers = 3.2) then
cylinder := plot3d( [2*cos(theta), 2*sin(theta), z], theta=0..2*Pi, 
z=0..6) :
plane := plot3d( (1-2*x +y), x=-2..2, y=-2..2) :
 with(plots) ;
display3d ( {cylinder, plane} ) ;  fi ; fi; end: