You can take the real and imaginary parts of the differential equation (in terms of the real and imaginary parts of the solution), and plot the results.  For example:

> de:= diff(z(t),t) = I*z(t) + t;
   realdes:= evalc({Re,Im}(eval(de,z(t)=x(t)+I*y(t))));
   Sol:= dsolve(realdes union {x(0)=0,y(0)=1}, numeric);
   plots[odeplot](Sol, [[t,x(t)],[t,y(t)]],t=0..4,colour=[blue,red]);

In other words, your constraint is

x[i,j] = 0 if x[u,v] = 0 for all u <> i, v > j

This kind of constraint is problematic for optimization because it leads to a feasible region that is not closed: the limit of a sequence of feasible points can be infeasible.  The result is that there will not exist
optimal solutions in which the constraint is binding.  In fact, suppose P is an optimal solution for this problem, but would not be an optimal solution without these constraints (C for short).  Then there is a solution P' that has a greater objective value that P, and is feasible for the problem without C, but does not satisfy C.  That would mean that in P', for some i,j we have x[i,j] > 0 but x[u,v] = 0 for all u <>i, v > j.  But if you increased all x[u,v] by a tiny bit, the constraints C would all be satisfied and the objective would still be greated than that for P, contradicting optimality of P.

Maple does not always display real solutions first, in fact in general it has no way of knowing which solutions are real.  In your example, the first solution happens to be real if k > 0, but there is no assumption that k > 0.

Perhaps something like this, suggested by <http://jap.aip.org/japiau/v32/i3/pS251_s1>:

> plot([sin(t+1/2),arctan(5*sin(t))/arctan(5),t=0..2*Pi]);

In a Plot component (rather than a regular plot) you can do this.  For example:

> de:= diff(y(x),x)=-y(x); ICs:= NULL;
P0:= DEtools[DEplot](de, y(x), x = -3 .. 3, y = -3 .. 3):

Now I insert a Plot component, right-click on it, click Component Properties.
Note that the Name is Plot0. 
Enter P0 for Plot Expression.
Check the box for "Make "execute code" the default manipulator".
Click on Edit... after "Action when Clicked", and change it to

global ICs;
local cx, cy;
use DocumentTools in
# Enter Maple commands to be executed when the specified
# action is carried out on the component.
# Use:
#    Do( %component_name );
# and
#    Do( %component_name = value );
# to set and get properties of the component.
# You can also use arbitrary expressions
# involving components, e.g.:
#    Do( %target = %input1 + 2*%input2 );
# Note the %-prefix to each component name.
# See ?CustomizingComponents for more information.
Do(cx = %Plot0(clickx));
Do(cy = %Plot0(clicky));
ICs := ICs, y(cx)=cy;
Do(%Plot0 = DEtools[DEplot](de, y(x), x = -3 .. 3, y = -3 .. 3, [ICs], linecolour=blue))
end use;

Click "OK", and again "OK".

Now when you click on the plot, a trajectory of the differential equation should appear, passing through the point you clicked on.

Rather than assigning a value to alpha, or unassigning it (both of which remove the previously assigned value that was the "assumed" version of alpha) you might try

>  assume(0 <= alpha);
>  f:= alpha^2:
>  eval(f, alpha = 2);

              4

You might also consider using "assuming" to make temporary assumptions on variables, rather than "assume" which makes a permanent change.

1) For example,

 > 2 + 3/4;
 

                                11/4

Maple automatically expresses the result as a common fraction, reduced to lowest terms.

2) For example:

> is( 4/7 > 9/23);

                 true

3) Since Maple automatically converts mixed fractions to common fractions, this can't be done directly.  However, you could do something like this:

> A:= 24/11:
   floor(A), A - floor(A);

                               2, 2/11

Actually plottools[arrow]([1,Pi/4], [1,0], 0.05, 0.1, 0.1) doesn't make the arrow go 1 unit to the right, it makes it go from [1,Pi/4] to [1,0], i.e. the [1,0] is interpreted as a point rather than a displacement.  If you wanted a displacement you should use a Vector <1,0> rather than [1,0]. 

You could use changecoords in the plots package to change coordinate systems.  However, note that

> plots[changecoords](plots[display](plottools[arrow]([1,Pi/4], <1,0>, 0.05, 0.1, 0.1)), polar);

is also not what you wanted: it makes the arrow go 1 unit in the radial direction, rather than 1 unit to the right (i.e. in effect the <1,0> is also interpreted in polar coordinates).  In addition, the arrow shape is somewhat strange: the width increases as it goes outwards.

The simplest work-around is to just convert the point from polar coordinates to rectangular coordinates.

There are really just 6 equations, not 9, because C . Transpose(C) is symmetric.  So you would expect at least 3 arbitrary variables

Note that if C is a solution, so is C U where U is any 3 x 3 orthogonal matrix (i.e. a matrix such that
U . Transpose(U) is the identity matrix).  The orthogonal matrices are a 3-dimensional manifold, so that accounts for your three arbitrary variables.  Moreover, the equation C . Transpose(C) = B . Transpose(B) (where these matrices are invertible) can be written as

B^(-1) . C . Transpose(B^(-1) . C) = B^(-1) . C . Transpose(C) . Transpose(B)^(-1)
   = B^(-1) . B . Transpose(B) . Transpose(B)^(-1) = IdentityMatrix

so B^(-1).C is orthogonal, i.e. any two solutions (in the invertible case) are related by an orthogonal matrix.

The Cholesky decomposition (in the case of a positive definite matrix Q, which yours is) finds a lower triangular matrix L such that L . Transpose(L) = Q.  Thus the general solution is C = L . U where U is an orthogonal matrix, and L can be obtained by

> L := LinearAlgebra:-LUDecomposition(Q, method=Cholesky);

                                    [2  0  0]
                                    [       ]
                               L := [2  2  0]
                                    [       ]
                                    [0  2  2]

By entering b:= a, you make both b and a point to the same object.  A change in b is then automatically a change in a.  If you want them to be different, try b := copy(a).

This problem is known as premature evaluation.  Your procedure func is designed to take a numerical argument  and return a numerical value.  But here you're calling it with a symbolic argument, and that doesn't work.  Besides the syntax plot(func, a .. b), a cure is to modify func so that it can accept a symbolic argument (in which case it returns unevaluated).  For example:

func := proc(T)
  if not type(T, numeric) then return 'procname(T)' end if;
 ... rest of code ....

In Unix you could use process[kill], but in most versions of Windows I don't think there is any way to do it.  This is a deficiency of the operating system.  There are utility programs available for this purpose, but I haven't had any experience with them.

A vector field A such that Curl(A) = B is called a vector potential of B.  You might try VectorPotential in the VectorCalculus package.

"not equal" in Maple is <>, not ~=.

Saving the worksheet does not save the values of variables.  You need to use the save command for that.  It does, however, save the results printed in the worksheet, and you may recover them using labels.  For example:

> A:= LinearAlgebra:-RandomMatrix(100,100);

The result, which looks like

                             [ 100 x 100 Matrix     ]
                             [ Data Type: anything  ]
                        A := [ Storage: rectangular ]
                             [ Order: Fortran_order ]

is given label (1).  If I save the worksheet (including the "references to large computation results") and re-open it, without executing any commands, the variable A is unassigned.  However, if I enter the label (1) followed by a semicolon and execute that, I get

                             [ 100 x 100 Matrix     ]
                             [ Data Type: anything  ]
                        A := [ Storage: rectangular ]
                             [ Order: Fortran_order ]

which restores the value of A.