> with(VectorCalculus):
  v := t -> <t^2, 3*t, t^3>;
  a:= D(v);

> a(t);
                                          2
                       2 t e  + 3 e  + 3 t  e
                            x      y         z

> a(5);

                         10 e  + 3 e  + 75 e
                             x      y       z

> BasisFormat(false):
   a(t), a(5);

                             [2 t ]  [10]
                             [    ]  [  ]
                             [ 3  ], [ 3]
                             [    ]  [  ]
                             [   2]  [75]
                             [3 t ]

Here it is using solve.

> Q:=solve({x1+x2+x3=1, x4+x5+x6=1, x7+x8+x9=1, x1+x4+x7=1, x2+x5+x8=1, x3+x6+x9=1});

Q := {x1 = -1+x5+x8+x6+x9, x2 = 1-x5-x8, x3 = 1-x6-x9, x4 = 1-x5-x6, x5 = x5, x6 = x6, x7 = 1-x8-x9, x8 = x8, x9 = x9}

Notice we have a parametric solution, with arbitrary values of 4 of the variables.  We iterate over all possible values
of those variables, checking whether the values of the other variables are in {0,1}, and print the resulting matrix in those cases.

> arb := convert(indets(map(rhs,Q)),list):
   iter:= combinat[cartprod]([{0,1}$4]):
   while not iter[finished] do
     X:= iter[nextvalue]();
     S:= map(t -> lhs(t) = subs(zip(`=`,arb,X),rhs(t)), Q);
     if map(rhs,S) subset {0,1} then print(subs(S,Matrix(3,3,[seq(x||i,i=1..9)]))) end if;
 end do:

Hints:

1) Draw a picture.  To see the shape of the curve, use more reasonable numbers than the ridiculous 224755712 and 369664.
2) Express the area as a function of the x coordinate of one vertex of the rectangle
3) Find the maximum of this function by the usual methods of calculus.

The trouble affects boundary value problems, but not the general solutions of differential equations.  So you can do this.

  for k from 1 to nmax+1 do
      gsol:= unapply(rhs(dsolve(des1[k], useint)), eta);
      F[k-1]:= subs(solve(eval({bcs(k-1)}, F[k-1]=gsol)), eval(gsol));
  end do;
      

It's not clear to me exactly what you're trying to do, but perhaps you mean something like this.

> U[0]:=  2*t - 3*t^2 + 3*x^2;
  for n from 1 to 4 do U[n]:= expand(U[0] + int(2*U[n-1]*diff(U[n-1],x),t)) end do;

Note that e is just another variable, not the constant  2.71828..., in Maple.  For that you want to use exp(1) (and e^2 is exp(2)).

Also, I think your differential equation in question 1 should be

y'' = (3 x y' - 4 y)/x^2

This particular one has no solution (e.g. because B is symmetric and full-rank, but A is not symmetric), but in general the way to solve it would be to write X with variables for entries, and use solve on the set of entries of X%^T . A . X - B.

See the help page ?numtheory,factorEQ.  The values of d for which the ring of integers of Z(sqrt(d)) is Euclidean are -1, -2, -3, -7, -11, 2, 3, 5, 6, 7, 11, 13, 17, 19, 21, 29, 33, 37, 41, 57, and 73.  -5 is not included.  IIRC, Z[sqrt(-5)] does not have unique factorization.

> evala(Norm(theta^2/2 - z));

-9*z-6*z^2-z^3+8

(the Maple tag is acting up again: that should be -9*z-6*z^2-z^3+8)

I suspect that all you need is the eval function, and maybe normal or simplify, but it's hard to say for sure without a specific example of what you're trying to do.

It's usually better to use gif rather than jpg for plots.

> plotsetup(jpeg, plotoutput="myfile.jpg", plotoptions="height=800,width=800");
   plot3d(x^2+y^2, x=-1..1, y=-1..1);
   plotsetup(default);

(or the same thing for .gif).  

Click on the plots for a closer look.

> ODEs:={4/5*F(x)*D(T)(x)-1/5*D(F)(x)*T(x)=(D@@2)(T)(x),  
 3/5*D(F)(x)^2-4/5*F(x)*(D@@2)(F)(x)=-(D@@3)(F)(x)+T(x)};
 BCs:= {F(0)=0,D(F)(0)=0,D(F)(10)=0,T(10)=0,D(T)(0)=-1};

Not surprisingly, dsolve fails to find a closed-form solution.  Try for a numeric solution.

> dsolve(ODEs union BCs, numeric);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

Following a suggestion in the help page ?dsolve,numeric,bvp, try continuation methods.  If the last boundary condition was D(T)(0)=0, there would be the trivial solution 0.

> Sol:=dsolve(ODEs union newBCs,numeric,continuation=lambda);

         Sol := proc(x_bvp) ... end proc

Plot this solution.

> plots[odeplot](Sol, [[x,T(x)],[x,F(x)]], x=0..10, colour=[red,blue]);

I don't know if there is a closed-form solution, but you might try a Laplace transform in one of the variables.

> with(inttrans):
   eval(laplace(PDE, x, s), {laplace(c(x,y),x,s) = U(y), c(0,y) = 39, D[2](c)(0,y) = 0});

s*diff(U(y),y)+A*s*U(y)-39*A+B*diff(U(y),y)+C*U(y)+D/s

This ODE in y can be solved explicitly.  The initial condition comes from the Laplace transform of BC2.  Note that you should always use exp(-x) rather than e^(-x).

> dsolve({%, U(0) = laplace(39*exp(-x),x,s)});

U(y) = exp(-(C+s*A)/(B+s)*y)*(39/(1+s)-(39*s*A-D)/(C+s*A)/s)+(39*s*A-D)/(C+s*A)/s

> Sol:= invlaplace(rhs(%),s,x);

Sol := 39*invlaplace(exp(-(C+s*A)/(B+s)*y)/(1+s),s,x)+D/C*invlaplace(exp(-(C+s*A)/(B+s)*y)/s,s,x)+(-A*invlaplace(exp(-(C+s*A)/(B+s)*y)/(C+s*A),s,x)*(39*C+D)+(39*cosh(1/2*C*x/A)*C-sinh(1/2*C*x/A)*(39*C+2*D))*exp(-1/2*C*x/A))/C

Unfortunately, in general Maple doesn't know the inverse Laplace transform of expressions involving exp(-(C+s*A)/(B+s)*y).  However, note that in the special case B = C/A, this evaluates nicely.

> S1:= eval(S, B = C/A);

S1 := 39/2*exp(-y*A-x)+1/2*D/C*exp(-y*A)+(-1/2*exp(-y*A-C*x/A)*(39*C+D)+(39*cosh(1/2*C*x/A)*C-sinh(1/2*C*x/A)*(39*C+2*D))*exp(-1/2*C*x/A))/C

which can be simplified a bit:

> map(simplify,collect(convert(S1,exp),exp));

1/2*D/C*exp(-y*A)+(39*C+D)/C*exp(-C*x/A)-1/C*D+39/2*exp(-y*A-x)-1/2*exp(-(A^2*y+C*x)/A)*(39*C+D)/C

Your fs is not the Laplace transform of ft: you have s/(s^2+4)-exp(-2*Pi*s)/(s^2+4), but it should be s*(1-exp(-2*Pi*s))/(s^2+4).

The Maple tag is acting up again: the first is s/(s^2+4)-exp(-2*Pi*s)/(s^2+4), the second is s*(1-exp(-2*Pi*s))/(s^2+4).

Maple can't fill 3-dimensional solid regions, but it can plot surfaces.  So you need to figure out the boundary of your region.  If you're plotting in a box, say,
-1 <= x,y,z <= 1, the boundary of the region x*z < y^2 will consist of the part of the surface x*z = y^2 in the box plus the parts of the planes x,y,z = (+/-) 1 for x^z < y^2.

> with(plots): with(plottools):
   eq := x*z-y^2;
  P0:= implicitplot3d(eq, x=-1..1, y=-1..1, z=-1..1, style=patchnogrid, axes=box, grid=[40,40,40]):
  display(subs(COLOUR(RGB,1,0,0)=NULL,[P0,
     seq(transform((x,y) -> [x,y,z])( 
       select(t -> (type(t,specfunc(anything,POLYGONS)) and has(t,COLOUR(RGB,1,0,0))),
          implicitplot(eq,x=-1..1,y=-1..1,filledregions=true))),z={-1,1}),
     seq(transform((x,z) -> [x,y,z])( 
       select(t -> (type(t,specfunc(anything,POLYGONS)) and has(t,COLOUR(RGB,1,0,0))),
          implicitplot(eq,x=-1..1,z=-1..1,filledregions=true))),y={-1,1}), 
     seq(transform((y,z) -> [x,y,z])( 
       select(t -> (type(t,specfunc(anything,POLYGONS)) and has(t,COLOUR(RGB,1,0,0))),
          implicitplot(eq,y=-1..1,z=-1..1,filledregions=true))),x={-1,1})]),
    shading=none,lightmodel=light3,orientation=[115,55]);