eq1 := -3.999168585*10^5*p1*q1^2-2.543619525*10^7*q1*p1^2+7.762255614*10^6*q1^3-3.999168584*10^5*p1^3-0.5000000000e-4*p1+0.1007754033e-3*q1 = 0:

eq2 := -3.999168585*10^5*q1*p1^2-7.762255614*10^6*p1^3+2.543619525*10^7*p1*q1^2-0.5000000000e-4*q1-3.999168584*10^5*q1^3-0.1007754033e-3*p1 = 0:

solve({eq1, eq2});

I did not understand that  pseudo rule means .

bcount:=n->nops(select(t->type(t, odd), [seq(binomial(n,k), k=0..n)])):

Example:

bcount(6);

                             4

It would be interesting to find an explicit formula for this.

Edited.

For any simple polygon (not necessarily convex) area can be found by shoelace formula . The algorithm implemented in Maple in the procedure  Area, which also solves  more general problem - symbolically (or numerically) finds the area of a figure, bounded by a non-selfintersecting piecewise smooth curve. 

convert(sin(x), FormalPowerSeries);

Should be

restart;

Ec := (1/2)*Kc*(2*Pi*h0/Pi/(a0^2+h0^2)-2*Pi*h/Pi/(a^2+h^2))^2;

simplify(Ec);

You have a system of 4 matrix equations in 4 unknowns matrices and matrix parameter y . I think that Maple can not solve the system symbolically, but Maple can solve it numerically for specific values of  B2 , B3  and  y (I took all the square matrix of the third order) . 

I have replaced the names of all the matrices with capital letters and converted the original system in 36 nonlinear scalar equations.

restart;

with(LinearAlgebra):

local D;

B2, B3, Y:=seq(RandomMatrix(3, generator=0..10), i=1..3);

A,B,C,D:=seq(Matrix(3, symbol=i), i=[a,b,c,d]); 

eq2 := A.B+C.D+A:

eq3 := A.C+C.D+C:

eq4 := A.B+C.A+B.C:

eq5 := A.B+A.D+B.C:

fsolve(map(t->Equate(op(t))[], [eq2=B2,eq3=B3,eq4=B2,eq5=Y]));

assign(%):

A,B,C,D;

Expr:=expand((x^2-x-3)^10);

select(t->sign(t)=-1, [op(Expr)]);

nops(%);

restart;

solve({-(1/3)*(eta+5)/(eta-3) = 3*eta/(-2+eta)});

assign(evalf(%[1]));

eta;

Instead of  LinearSolve  should be  LinearAlgebra[LinearSolve] 

If you have  m=1..2  and  x=0..2  then it is a function of 2 variables, so

plot3d(sin(m*x), m = 1 .. 2, x = 0 .. 2, axes=normal);

restart;

taylor(f(x), x = gamma);

subs([x-gamma=e[n], seq((D@@k)(f)(gamma)=k!*c[k]*f(gamma), k=1..5)], %);

 Addition. May be should be  f[k](gamma)  instead of  f(gamma) ?

In addition to the plotting of the ellipse was asked to find a canonical equation. The simplest way to do this - using of the formulas from here

restart;

local D;

f:=(x,y)->3*x^2-3*x*y+6*y^2-6*x+7*y-9;

coeffs(f(x,y));

A,B,C,D,E,F:=%;   # Assigning names to coefficients

theta:=1/2*arctan(B/(A-C));   # Finding the angle of rotation of the canonical system of coordinates relative to the original  coordinate system

solve({-2*A*xc-B*yc=D, -B*xc-2*C*yc=E});   # Finding of the center of the ellipse

assign(%);   # Assigning names to coordinates of the center

# The connection between the canonical  coordinates  and the original  coordinates

x:=xcan*cos(theta)-ycan*sin(theta)+xc;   

y:=xcan*sin(theta)+ycan*cos(theta)+yc;

Eq:=simplify(expand(f(x,y)));   # The equation of the ellipse in canonical coordinates

xcan^2/simplify(sqrt(1/coeff(Eq,xcan^2)*(-tcoeff(Eq))))^`2`+ycan^2/simplify(sqrt(1/coeff(Eq,ycan^2)*(-tcoeff(Eq))))^`2`=1;  # Canonical equation

Visualization:

xcan:=plot(yc+tan(theta)*('x'-xc), 'x'=-2..3.5, color=black):

ycan:=plot(yc-1/tan(theta)*('x'-xc), 'x'=0.1..1.5, color=black):

Ellipse:=plots[implicitplot](f('x','y'),'x'=-2..3.5,'y'=-2..1.5,color=red,thickness=2, gridrefine=5):

labels:=plots[textplot]([[0.4,1.3,"ycan"],[3.2,0.75,"xcan"]], font=[TIMES,ROMAN,14]):

plots[display](xcan,ycan,Ellipse,labels,scaling=constrained);

plot(f(t,a0,b0,c0),  t=0..10);  # or

plot(eval(f(t,a,b,c), {a=a0, b=b0,c =c0), t=0..10);

restart;

Puis:=proc(X, n::integer)

option remember;

if n=0 then 1 elif

n>0 then Puis(X, n):=Puis(X, n-1)*X  elif

n<0 then Puis(X, n):=Puis(X, n+1)/X  fi;

end proc:

Examples of use:

Puis(4,0),  Puis(4,3),  Puis(4,-3),  Puis(4,6),  Puis(4,-6);

                           1, 64, 1/64, 4096, 1/4096

Addition: I understood the meaning of the problem is not to use the sign  ^

rhs  is superfluous. Should be

a := .1994;

modfit3 := a*x^1.5;

f := unapply(modfit3, x);

Compare:

rhs(b);

   Error, invalid input: rhs received b, which is not valid for its 1st argument, expr