Rather than using that package, you can use orthopoly[P] for the Legendre polynomials.  This should be quite fast.  For example:

> Leg:= orthopoly[P];
   In:= (f,g) -> int(f*g,x=-1..1);

# check the normalization

> seq(In(Leg(n,x),Leg(n,x)), n=0..5);

      2 = 2, 2/3 = 2/3, 2/5 = 2/5, 2/7 = 2/7, 2/9 = 2/9, 2/11 = 2/11

> LegSeries:= (f,n) -> add(In(f,Leg(j,x))*Leg(j,x)*(j+1/2),j=0..n)  ;

> F:= piecewise(x<0,0,1);
   plot([F, LegSeries(F,20)], x=-1..1);

You can also prove uniqueness using a single objective: the sum of the variables (including slack variables) that are 0 in the initial solution.  With NewC as in my previous posting:

> nbz:= map(t -> (t=0), indets(map(rhs,NewC)));
  

> xzeros:= select(t -> (eval(eval(t,NewC),nbz)=0),{seq(cat(x,i),i=1..7)});
   azeros:= select(t -> (eval(eval(t,NewC),nbz)=0), {seq(eval(cat(a,i)),i=1..7)});

>  obj:= convert(xzeros,`+`) - convert(azeros,`+`); 

> simplex[maximize](obj,cons,NONNEGATIVE);

> eval(obj,%);

    0

You have a system of linear inequalities, which may be considered as the constraints of a linear programming problem.  If you just want one numerical solution, you might try something like this:

> cons:= {b = 1, a1 <= 0, a2 <= 0, a3 <= 0, a4 <= 0, a5 <= 0, a6 <= 0, a7 <= 0};
   Optimization[Maximize](0,cons,assume=nonnegative);

[0., [x1 = 0., x2 = 0., x3 = .625000000000000e-1, x4 = .312500000000000, x5 = .250000000000000, x6 = .312500000000000, x7 = .625000000000000e-1]]

In fact this turns out to be the only solution, but that's not obvious here.

Or you could try the simplex package.

> simplex[feasible](cons,NONNEGATIVE,'NewC');

                           true

> NewC;

Here _AR is an artificial variable, which must be 0 for feasibility; _SL1, _SL2 etc are slack variables that must be >= 0.  Notice that the equation for x2

can only be satisfied if all the variables it contains are 0. Setting them to 0:

> subs(map(t -> (t=0), indets(NewC[3])), NewC);

[_SL2 = 13/16, x7 = 1/16, 0 = 0, x3 = 1/16, x4 = 5/16, x5 = 1/4, _SL8 = 1, x6 = 5/16]

Note that the rather complicated solution returned by pdsolve (which is not the general solution) only depends on the variables x,y,z through _C1*x + _C2*y + _C3*z.  Here _C1, _C2 and _C3 are arbitrary constants.  The pde is invariant under rotations, so we may as well look for solutions where f is a function of x. 

> dsolve(eval(pde, f(x,y,z) = f(x)));

,

That's what I did to get rid of the square root.  And the result is

Your differential equation has symbolic parameters a,b,c,d,e as well as an additional unknown function z(x).  You also didn't specify initial or boundary conditions.  Without knowing what these are (or at least most of them), it's going to be difficult to determine the maximum.

Yes, the normal should be perpendicular to the tangent, and that's exactly what your graph shows: the tangent in green, the normal in gold.  It looks fine on my computer.  I don't know what's bothering you about this.

To make sure the scales are the same on both axes, you could use the scaling=constrained option in your plot command.

I think what you mean is f(x) = (6*x^2+5*x-4)/(2*x^2-7*x-4). Since you used Tangent(f(x),x=-3), you probably want the tangent at (-3,1), not (0,1).
And for the normal, you mean 35/12*x + 39/4.

If you have doubts about the answer, check it yourself: is the slope correct?  Is the value at x=-3 correct?

What do you mean by the cross product of two matrices? Cross product is an operation on vectors in R^3, not on matrices.

Like this?

N:= 10; 
S[1]:= -gamma; 
for nn from 2 to N do S[nn]:= S[nn-1]+1/(nn-1); end do;
step:= n -> ([n,S[n]],[n+1,S[n]]);
theHarmonicNumbersPlot:= plot(step~([$2..N-1])):
theHarmonicFunctionPlot:= plot(Psi(n), n=2..N, colour=black):
plots:-display([theHarmonicFunctionPlot,theHarmonicNumbersPlot]);

You have 3 variables x,y,z, and four inequalities, determining a 3-simplex in R^3, not a 2-simplex in R^2.
If you replaced one inequality by an equality, you'd have a 2-simplex, but still in R^3.  In R^2, the feasible region for a system of inequalities can be plotted using inequal in the plots package. 

I might do the dynamic programming this way (assuming the problem is to select a subset of [$1..N] with total weight exactly W and maximum value).  Vmax[n,w] is the maximum value for the first n items with total weight w, and Smax[n,w] is true if one of those optimal solutions includes item n.  At the end, Q is the optimal selection from the N items for total weight W.



> Vmax[0,0]:= 0:
 for y from 1 to W do
       Vmax[0,y]:= -infinity;
 end do:
 for n from 1 to N do
      for y from 0 to w[n]-1 do
         Vmax[n,y]:= Vmax[n-1,y]; Smax[n,y]:= false     
      end do;
      for y from w[n] to W do
         v1:= Vmax[n-1,y];
         v2:= v[n] + Vmax[n-1,y-w[n]];
         if v1 > v2 then Vmax[n,y]:= v1; Smax[n,y]:= false
         else Vmax[n,y]:= v2; Smax[n,y]:= true
         end if
      end do
   end do:
   Q:= NULL; y:= W;
   for n from N to 1 by -1 do
     if Smax[n,y] then Q:= Q,n; y:= y-w[n]; end if;
   end do:
   Q:= {Q}:
   `Optimal value`, Vmax[N,W], `for total weight`,W,`with selection`, Q;

Walking the structure once is not too bad either.

>C:= rtable_scanblock(M,[1..20,2..2],
      proc(v,i,r) if v < r[2] then r elif v = r[2] then [[r[1] ,i[1]],v]
         else [[i[1]],v] fi end proc,
      [{1},M[1,2]], 
      proc(v,i,r) if v > r[2] then r elif v = r[2] then [[r[1],i[1]],v]
         else [[i[1]],v] fi end proc,  
      [{1},M[1,2]]): 
 A:= indets(C[1][1],integer);
 B:= indets(C[2][1],integer);

I would avoid the geometry and geom3d packages, which I find rather cumbersome.  For drawing triangles, it's simplest to just use polygonplot in the plots package.  For example:

> with(plots): polygonplot([[1,2],[3,4],[0,0]], colour=red, scaling=constrained);

For tetrahedra, use polygonplot3d for each face.  For example:

> Vertices:= [[1,2,0],[2,3,1],[3,2,1],[4,2,3]];
   display([seq(polygonplot3d(subsop(i=NULL,Vertices)),i=1..4)],scaling=constrained);

For (a):

PyramidProj:= proc(a,h,A,B,C)
    uses LinearAlgebra;
    local V,verts;
    V:= GramSchmidt([<A,B,C>,<1,0,0>,<0,1,0>,<0,0,1>], conjugate=false,normalized=true);
    verts:= map(v -> [DotProduct(v,V[2],conjugate=false),
                               DotProduct(v,V[3],conjugate=false)],
       [<0,0,0>,<a,0,0>,<a,a,0>,<0,a,0>,<a/2,a/2,h>]);
    simplex[convexhull](evalf(verts),output=area);
  end proc;

For an approximation of (b), you might try

ConeProj:= proc(r,h,A,B,C)
    uses LinearAlgebra;
    local V,verts;
    V:= GramSchmidt([<A,B,C>,<1,0,0>,<0,1,0>,<0,0,1>], conjugate=false,normalized=true);
    verts:= map(v -> [DotProduct(v,V[2],conjugate=false),
                               DotProduct(v,V[3],conjugate=false)],
       evalf([<0,0,h>,seq(<r*cos(2*Pi*t),r*sin(2*Pi*t),0>,t=0 .. 0.999,.001)]));
    simplex[convexhull](verts, output=area);
  end proc;