One way is to extract the data from the plot.

sol1:= dsolve(
     [diff(u(theta),theta) = 427.2461*u(theta)+385620.123/u(theta)-25671.3871,
      u(0) = .6
     ], numeric
):
plots[odeplot](sol1, 0..0.18, color = red):
A:= op([1,1], %);

A[1..10, 1..2];

There are two ways to "clear" a variable.

1.

unassign('a');

2.

a:= 'a';

In either case the single quotes are required.

@verdin The solution for the real and for the imaginary parts could be run in separate processes. This has the potential to cut your time in half. There is no obvious way to parallelize the individual solves, so I don't think that you could get any more savings than that.

Make your two equations to be solved into a list and apply Grid:-Map to it. This will run the solves in separate processes.

How about this. The key to emphasizing the inequality aspect is using plot option filled. I overlayed two copies of the plot so that I could have both the grid and the surface without the grid being black. The transparency option also serves to highlight the "solid" nature of the object being plotted (i.e., that it is more than a surface).

Fplot:= x^2*y, x= 0..1, y= 0..1, grid= [10$2], filled:
plots:-display(
     [plot3d(Fplot, style= wireframe), plot3d(Fplot, style= patchnogrid)],
     orientation= [130,65], axes= boxed,
     view= [(0..1) $ 3], transparency= 0.2
);

The mapping < 2, 2 > -> < 1, 1, 1 > and < 1, 4 > -> < 1, 2, 0 > is just a random choice made to make an example. The point of the wiki page is that once those two choices are made, the action of the mapping on any other vector in R^2 is determined.

It's a trivial application of LinearAlgebra:-LinearSolve:

d:= < <1,0,0> | <0,-2,0> | <1,0,1> >: #The basis from the wiki page.
Rep:= (D,v)-> LinearAlgebra:-LinearSolve(D,v);
Rep(d, <1,1,1>);

Here's an example to show what the results returned by DirectSearch:-SolveEquations mean and how to extract the individual results.

Download SolveEquations.mw

You need to split the inequality into two and use solve. Like this:

T:= 4477.25+.25*(t-32450):
solve({4477.25 <= T,T <= 16042.25}, t);

Use * to represent multplication:

y - y1 = m*(x - x1);

Your original formulation returns the error "Initial Newton iteration is not converging":

restart:
ode:=
     diff(f(eta),eta$3) + f(eta)*diff(f(eta),eta$2) - (diff(f(eta),eta))^2 +
     lambda*(2*f(eta)*(diff(f(eta),eta))*(diff(f(eta),eta,eta)) -
            f(eta)^2*diff(f(eta),eta,eta,eta)) -
     k*(diff(f(eta),eta)-1) =
     0:
params:= [lambda= 0.1, k= 0.2, inf= 10]:
BCs:= f(0) = 0, D(f)(0) = 1, D(f)(inf) = 0:
dsolve(eval({ode, BCs}, params), numeric);

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

This is a common error for numeric BVPs. It can often be corrected by introducing a "continuation parameter." (See ?dsolve,bvp,advanced .) Notice the C that I added to your ode below. It is in front of the fourth term, next to lambda. There is a lot of guesswork involved in deciding where to put the parameter. I usually start with making it a coefficient of the most complicated term, which is what worked here.

ode:=
     diff(f(eta),eta$3) + f(eta)*diff(f(eta),eta$2) - (diff(f(eta),eta))^2 +
     C*lambda*(2*f(eta)*(diff(f(eta),eta))*(diff(f(eta),eta,eta)) -
            f(eta)^2*diff(f(eta),eta,eta,eta)) -
     k*(diff(f(eta),eta)-1) =
     0:

Sol:= dsolve(eval({ode, BCs}, params), numeric, continuation= C);
plots:-odeplot(Sol, [eta, f(eta)], 0..10, thickness= 2);

I is reserved as the imaginary unit in Maple. You can change this by declaring another symbol to be the imaginary unit, let's say II, by the command

interface(imaginaryunit= II);

In Maple 17+, you can also declare local variables at the top level:

local I;

and then use I as you wish.

The imaginary parts are small and probably spurious (i.e., due simply to round-off error). You can plot the real parts with Re:

plot(Re(eq), x= 0..1);

If you have some theoretical reason to believe that this expression eq should be purely real for x= 0..1, then you can verify that the imaginary parts are truly small and safe to ignore by plotting the imaginary parts:

plot(Im(eq), x= 0..1);

Doing this, I see that the real parts are of order 1 and the imaginary parts are of order 10^(-6). Considering the huge size of your coefficients, I think that it is safe to ignore the imaginary parts, if you have some theoretical reason to do so. 

Student:-Calculus1:-VolumeOfRevolution(
     x^2, sqrt(x), x= 0..1, axis= horizontal , distancefromaxis= -1,
     output= animation, partition= 40,
     volumeoptions= [color= yellow]
);

If you you want to do the rotation about the line x = -1 by the washer method, then you need to integrate with respect to y. Because your functions are inverses of each other, the command is the same as above with x changed to y, which in reality changes nothing. More interesting is to do it by shells. For this change the axis option to axis= vertical.

To get a Maple plot similar to the kgraph plot, use the view option to restrict the view:

plot(A[.., 1], A[.., 2], view= [-1.5e-9..2.5e9, -2..2]);

Since you are using (1D) Maple Input, it is incorrect to define a function (procedure) by

g(t):= f(x(t), alpha(t), beta(t))

The correct syntax is

g:= t-> f(x(t), alpha(t), beta(t)).

If you use the "g(t):=", then you are defining an expression that will only work when the t is explicitly used. That is, g(t) will work, but g(1.0) will not.