If you want the solution with F > 0, you might be able to get that by specifying the inequality F > 0 together with the equation to solve.  Thus:

> solve({ K = S*F^n, F > 0}, F);

Or if you want the real solutions, try using the inequality F > -infinity.

Perhaps something like this?

> with(Statistics):
   X:= Sample(Uniform(0,1),30);
   Y:= Sample(Normal(0,1), 30);
   Data := zip(`[]`,X,Y);

This requires constructing the PLOT data structure yourself:

Data:= op(fscanf("8441_pvtszfrmt.txt","%{41,91,3}fa"));
 zmax,zmin:= (max,min)(seq(seq(Data[i,j,3],j=1..91),i=1..41));
 deg:= evalf(Pi/180);
 PLOT(MESH([seq([seq([Data[i,j,1]*cos(Data[i,j,2]*deg),
    Data[i,j,1]*sin(Data[i,j,2]*deg)],j=1..91)],i=1..41)],
 COLOR(HUE,seq(seq((Data[i,j,3]-zmin)/(zmax-zmin),
    j=1..91),i=1..41))),
 STYLE(PATCHNOGRID), AXESSTYLE(BOX), SCALING(CONSTRAINED));

Your differential equation is a first-order linear ODE, and the standard solution formula applies:
 

> dsolve(K/i(t)*C*diff(v(t),t)+v(t)=i(t), v(t));

v(t) = (Int(i(t)^2*exp(1/K/C*Int(i(t),t))/K/C,t)+_C1)*exp(Int(-i(t)/K/C,t))

Of course, there's no guarantee that the integrals can be done in closed form, and indeed in the case i(t)=sin(w*t) you get an integral that probably does not have a closed form:

v(t) = exp(1/K/C/w*cos(w*t))*int(-1/2*(-1+cos(2*w*t))*exp(-1/K/C/w*cos(w*t))/K/C,t)+exp(1/K/C/w*cos(w*t))*_C1

As for "heavyside", in Maple that function is called Heaviside.

It seems the <maple> tag has stopped working again 

Your f is missing a ")".  I assume it should be

f := (u, v) -> (u*(1-u))-a*(u*v/(v+u));

But the main problem is with g, where you are missing multiplication signs after
the first v and the a: this causes Maple to treat v and a as functions.
Using

g := (u, v) -> b*v*(a*(u/(v+u))-d);

you get the steady states

{v = 0, u = 1}, {v = -(-1+a-d)*(a-d)/d, u = 1-a+d}

I suspect you want

> dsolve({deq2, ci}, x(t));

not

> dsolve({deq2,x(t)});

(which actually tells dsolve that you want one of the equations to be x(t) = 0).

I assume your equations involve polynomials.

This might be overkill, but:

In the Groebner package, use IsZeroDimensional to check whether the set of equations has a finite set of solutions.  Use IsProper to check whether there is at least one solution (over the complex numbers).  Use MaximalIndependentSet to find a maximal set of algebraically independent variables.  In your example:
 

> polys := {2*x+3*y, y-2+z};
   with(Groebner):
   IsZeroDimensional(polys);

                      false

> MaximalIndependentSet(polys);

                    { y }

By the way, what "we all know" is that n equations are necessary, but are not necessarily sufficient.

It stops with an error because you're trying to assign a value to A, which is a formal parameter of the procedure sub1.  That's a bad idea.
I think what you want is

sub1:= proc(i, k, f, A)
      local L;
        if   i <= f then  A[i]:=A[i]+1; L:=convert(A, `list`);
                 if sum(L[j]*q^j, j=1..f)<=  k then
                     print(L);  sub1(1, k, f, Array(L));
                  else A[i]:=0; sub1(i+1, k, f, A);
                 fi;
        else print("DONE");
        fi;
   end;

Then lis(13) produces the  output

                              [1, 0, 0]
                              [2, 0, 0]
                              [3, 0, 0]
                              [4, 0, 0]
                              [5, 0, 0]
                              [6, 0, 0]
                              [0, 1, 0]
                              [1, 1, 0]
                              [2, 1, 0]
                              [3, 1, 0]
                              [4, 1, 0]
                              [0, 2, 0]
                              [1, 2, 0]
                              [2, 2, 0]
                              [0, 3, 0]
                              [0, 0, 1]
                              [1, 0, 1]
                              [2, 0, 1]
                              [0, 1, 1]
                                "DONE"

  But you reallly don't need to keep converting between lists and Arrays.  A simpler version is

sub1:= proc(i, k, f, A)
       if   i <= f then  A[i]:=A[i]+1; 
                 if add(A[j]*q^j, j=1..f)<=  k then
                    printf("[%{c,}d]\n",A); sub1(1, k, f, A);
                 else A[i]:=0; sub1(i+1, k, f, A);
                fi;
       else printf("DONE");
       fi;
  end;

Use fprintf.  For example:

> A := << a^2| 3.5>, < 2 | b>>;
  fprintf("myfile.txt","%{c,}a\n", A);     
  close("myfile.txt");

and the file should contain:

a^2,3.5
2,b
 

Actually I think surfdata is the command you want.  Note that I'm converting the angle measurements (which seem to be in degrees) to radians.

> Data:= op(fscanf("8441_pvtszfrmt.txt","%{41,91,3}fa"));
   for i from 1 to 41 do 
      for j from 1 to 91 do
          Data[i,j,2]:= Data[i,j,2]*Pi/180
    end do
    end do:
   with(plots):
   surfdata(Data, coords=cylindrical, axes=box,     
     style=patchcontour,
     orientation=[0,180], shading = zhue);
  

First of all, it's best not to use I for that, because I in Maple is sqrt(-1).

The usual Maple command for Laplace transform is laplace in the inttrans package, but that does not help here.  I suspect there may be no closed form
for this Laplace transform.  You could of course use the definition of the Laplace transform as an integral:

> L := int((1-exp(-sqrt(sin(w*t))))/sqrt(sin(w*t))*exp(-s*t),t = 0 .. infinity);

Again, no closed form is found.  Numerical evaluation is possible in principle but will be slow because of all the singularities.  However, we can try for a version that will be easier to evaluate numerically.

First of all, nondimensionalize: L = L1/w where L1 depends only on s/w.  So in effect we may assume w=1.

> L1 := eval(L, w=1);

Now the integration can be restricted to the interval [0, Pi/2] by the symmetries of sin(t).  And then the substitution t = r^2 gets rid of the singularity at t=0.  The end result is this, which evaluates relatively quickly for small s:

> L2 := proc(s) 
 1/(1 - exp(-2*Pi*s))*(
 Int(2*r*(exp(-s*(r^2+Pi))*sin(sin(r^2)^(1/2))+exp(-s*r^2)
    -exp(-s*r^2-sin(r^2)^(1/2))+exp(-2*Pi*s+s*r^2)
    *sin(sin(r^2)^(1/2))+exp(-s*(Pi-r^2))-exp(-sin(r^2)^(1/2)
    +s*(-Pi+r^2)))/sin(r^2)^(1/2),
   r = 0 .. sqrt(1/2*Pi))
 +I*Int(2*r*(-1+cos(sin(r^2)^(1/2)))*(exp(-s*(r^2+Pi))
    +exp(-s*(2*Pi-r^2)))/sin(r^2)^(1/2),
   r = 0 .. sqrt(1/2*Pi))) 
 end proc;

arctan(b/a) is not the same as arg(a+b*I).  For example, consider a=b=-1.

On the other hand, you could use the two-variable form of arctan:
arg(a+b*I) = arctan(b,a).

Maple seems to be using the method of variation of parameters to get the particular solution of this inhomogeneous linear ODE.   With basic solutions y_1(x) = sin(2 x) and y_2(x) = cos(2 x) of the homogeneous equation, the standard VOP formula gives the particular solution 1/8*cos(2*x) + x/4*sin(2*x)

Hints:

1) don't forget * for multiplication (or space if you're using 2D math input)

2) diff(expr, x$n) is the n'th partial derivative of expr with respect to x.

3) sum produces a summation (finite or infinite)

4) However, your sum doesn't converge

The main problem, I think, is that t is a global variable, and the t's in the different functions interfere with each other.  Here's a simpler example.

> Q:= proc(f)  x -> int(f(t)*(t-x), t = x .. x+1) end proc;

> Q(f)(s);

int(f(t)*(t-s),t = s .. s+1)

But:

> Q(f)(t);

0

Why?  Because when you substitute t for x in t-x you get t-t  = 0. 

And therefore

> Q(Q(f))(x);

0

because Q(Q(f)) calls Q(f)(t).

Of course you want the t in the definition of Q to be a "dummy variable", not the same as the t that you substitute for x.  But that doesn't happen automatically: the only way to make it so is to make t a local variable.  Thus:

> Q:= proc(f) local t;  x -> int(f(t)*(t-x), t = x .. x+1) 
      end proc;

> Q(f)(t);

int(f(t)*(t-t),t = t .. t+1)

You can't tell by looking at it which t is which, but Maple knows...

 
> Q(Q(f))(x);

int(int(f(t)*(t-t),t = t .. t+1)*(t-x),t = x .. x+1)