You can rewrite your equation as a system dx/dt = f(x,y), dy/dt=g(x,y) e.g. as follows:

restart;
with(DEtools):
a:=.4384477959e16; b:=.133e17; c:=.110224e16; d:=.219223898e16;
#Scaling the constants
S:=1e-15:
sys:=expand({diff(y(t),t)=S*sqrt(y(t))*(a*sqrt(x(t))-b*y(t)^(1/4)),diff(x(t),t)=S*x(t)^2*(c*y(t)^(7/4)-d*sqrt(x(t)))});
DEplot(sys, [x(t),y(t)], t = -1 .. 1,x=0..5, y = 0 .. 5, [[x(0)=1,y(0) = 2], [x(0)=2,y(0) = 2], [x(0)=3,y(0) = 2],[x(0)=4,y(0) = 2]], linecolor = [green, blue, yellow,orange], arrows = medium, size = 1, stepsize = 0.01 );

Fe := p->fsolve(1-x+p = x, x);
plot(Fe,0..1);
plot('Fe(a)',a=0..1);

The reason for the error message from plot(Fe(a),a=0..1) (without quotes) is that arguments to plot are evaluated before plot sets in. Thus Fe(a) is evaluated, but fsolve(1-x+a=x,x) cannot be solved unless a is a number.

Fe evaluates to the name Fe, 'Fe(a)' evaluates to Fe(a), thus no problem is perceived at that stage in those cases.

In some situations you may want to define your function so that it returns unevaluated for non-numeric input, e.g. like this:

F2:=p->if not type(p,realcons) then 'procname(p)' else fsolve(1-x+p = x, x) end if;
plot(F2(a),a=0..1);

Here F2(a); just returns F2(a) (if a is just the name a).

Absolutely nothing is imposed on x by the assignment eq := x = a+x^b;

You could try assume(eq); but since that involves a and b on which you already made assumptions you need to use additionally.

First consider a modified and much simpler problem.

restart;
assume(a > 0);
assume(b < 1);

eq := x = a+x*b;
additionally(eq);
about(a);
about(x);
is(x>0);

Now your original problem, where you get the answer FAIL from is(x>0) and coulditbe(x>0).

restart;
assume(a > 1);
assume(b > 1);

eq := x = a+x^b;

additionally(eq);
about(a);
about(x);
is(x>0);
coulditbe(x>0);


Answering the question by zhong chen:

When using piecewise I don't see any possibility for different styles for the component expressions.

But you could do like this:

f1:=piecewise(x < 0, sin(x),undefined);
f2:=piecewise(x < 0,undefined, x<2*Pi,cos(x),undefined);
f3:=piecewise(x < 2*Pi,undefined,2);
plot([f1,f2,f3],x=-5..10,color=[red,blue,black],thickness=3,linestyle=[1,2,3]);

Inspired by an example in the help page for overload here is a start.

Suppose you want the union of two lists to return the list obtained by appending the second one to the first.

Union := module() option package;  export `union`;
    `union` := proc(a::list, b::list) option overload;
        [op(a),op(b)];
    end proc;
end module;

with(Union);
{7,8} union {97,7,8,9};
[7,8] union [97,7,8,9];
[7,8] union [97,7,8,9] union [a,b];
#The following is not handled and results in an error:
{7,8} union [97,7,8,9];

I wouldn't try adding the results of two DEplot graphs. The numeric procedure adjusts the stepsize so you don't necessarily have the same number of points in the two graphs.

A much better idea is to use dsolve/numeric and odeplot as in the following example, where the odes are coupled, but that is irrelevant.

with(plots):

ode1:=diff(x(t),t)=sin(x(t))+cos(y(t));
ode2:=diff(y(t),t)=sin(x(t)*y(t));
sol:=dsolve({ode1,ode2,x(0)=1,y(0)=1},numeric);
odeplot(sol,[t,x(t)],-10..10);
odeplot(sol,[t,y(t)],-10..10);
odeplot(sol,[t,x(t)+y(t)],-10..10);
#plot in phase plane:
odeplot(sol,[x(t),y(t)],-10..10);

However, if you insist, and you check that the list of t-values are the same for the two plots, then you could do as follows.

with(DEtools):
DEplot({ode1,ode2},[x(t),y(t)],t=-10..10,[[x(0)=1,y(0)=1]]);
DEplot({ode1,ode2},[x(t),y(t)],t=-10..10,[[x(0)=1,y(0)=1]],scene=[t,x(t)]); p1:=%:
DEplot({ode1,ode2},[x(t),y(t)],t=-10..10,[[x(0)=1,y(0)=1]],scene=[t,y(t)]); p2:=%:
L1:=op(indets(p1,listlist));
L2:=op(indets(p2,listlist));
nops(L1);
nops(L2);
map2(op,1,L1-L2);
L:=zip((u,v)->[u[1],u[2]+v[2]],L1,L2);
plot(L);

latex(Matrix(m1),"C:/users/desktop/work/m1.tex");

Notice that you should use either \\ or /.

The answer is probably not so surprising if you try the following, where all constants are set to 1:

solve( eval({eqn1,eqn2,eqn3},{M=1,a=1,b=1,c=1,d=1}),{x,y,z} );

evalf(%);
     {x = 0.9476095982, y = 0.9476095982, z = 0.9476095982}

You do get an answer from solve, but it involves RootOf's which probably cannot be expressed in terms of known functions.

For given values of the constants you can solve the system numerically using fsolve or (as done above) by solve followed by evalf.

p:=plot([sin(x), cos(x)],x=0..1, legend=["sin","cos"]):
L:=indets(p,listlist):
nops(L);
type(L[1],listlist);
indets(p,specfunc(anything,LEGEND));

Exact solution on a less trivial initial value problem is unlikely to result in an answer.

However, numerical solution works fine on the original system.

(The crucial problem in your worksheet was the presence of ,bigsys at the end of the dsolve command).

restart;
with(plots):
#Changing the initial value for x to make the solution less trivial:
res:=dsolve({x(0)=1, y(0)=0, z(0)=-1, diff(x(t), t)=I*x(t)+y(t), diff(y(t), t)=-I*y(t)+x(t)*z(t), diff(z(t),t) = -2*(conjugate(x(t))*y(t)+x(t)*conjugate(y(t)))},numeric,output=listprocedure);
X,Y,Z:=op(subs(res,[x(t),y(t),z(t)]));
X(1.2345);
odeplot(res,[t,Re(x(t))],0..5);
complexplot(X(t),t= 0..5,axes=boxed,thickness=2,color=black,legend=[x(t)]);
complexplot(Y(t),t= 0..5,axes=boxed,thickness=2,color=black,legend=y(t));
#If z(0) is real, z(t) will be real for all times t because of the ode for z:
plot(Z(t),t= 0..5,axes=boxed,thickness=2,color=black,legend=Z(t));

Executing your worksheet I get the same output as you do.

However, initially when I just copied from your post and pasted it into Maple I got what I expected.

After that using F3 I tried separating your worksheet into 2 execution groups, the first one ending with alpha1:= 1-w1-w2:

and the second one starting with interface(...).

Then your worksheet behaved as expected.

The problem could be that interface() has to do with just that, the interface, so that execution separation is not irrelevant.

But after uniting the two execution groups again using F4, I still get the desired output, so actually I don't know what is going on.

To make Maple write a particular expression in another particular way can be difficult.

If you already know the answer it is often a waste of time, but it  does present some challenges, if that is what you are looking for.

Easier is it to simplify the difference.

In this case you could do like this:

res1:=convert(sin(x)^4/cos(x)^2,tan);
simplify(`*`(op(1..3,res1)))*op(4,res1);

The last step brings to mind Robert Israel's procedure applysome, which allows for applying a function/procedure to some subset of the operands of a sum or  product.

Another way:

convert(sin(x)^4/cos(x)^2,csccot);
eval(%,{csc=1/sin,cot=1/tan});

Change the lines defining nAg and kAg to

nAg := unapply(Spline(nAgData, lambda), lambda);
kAg := unapply(Spline(kAgData, lambda), lambda);

Then Spline(nAgData, lambda) and Spline(kAgData, lambda) are first evaluated, then made into functions of the variable (i.e. name) lambda.

with(plots):

spacecurve(<cos(x), sin(x), exp(sin(x))>,x=0..2*Pi,axes=boxed,thickness=4):
spacecurve(diff~(<cos(x), sin(x), exp(sin(x))>,x),x=0..2*Pi,axes=boxed):
spacecurve(int~(<cos(s), sin(s), exp(sin(s))>,s=0..x),x=0..2*Pi,axes=boxed,thickness=2):
display(%%%,%%,%);

p:=proc(fx::algebraic,rx::name=range(realcons))  local x,c;
   x:=lhs(rx);
   c:=diff(fx,x,x)/(1+diff(fx,x)^2)^(3/2);
   plot([fx,c],rx,_rest)
end proc;
p(x^2,x=-2..2);
p(cos(x),x=-Pi..Pi,color=[black,blue]);
p(sqrt(1-x^2),x=-1..1,thickness=2);

The procedure accepts additional arguments besides the expression and range. They are passed on to plot in the form of _rest .

Comment: I just noticed that this question was posed twice and already answered in its first version.

Any reason to ask the same question twice?