It works for me:

> DocumentTools[SetProperty](ComboBox0,itemList,[x1, x2, x3, x4, x5, etc, etc]);

I have some doubts that your program does what you intend, but you might start with this.  Note that the coefficient of p[j]*(x-I)^k in the Taylor series of f*Q-P about x=I is the coefficient of (x-I)^k in the Taylor series of f*x^j about x=I.  Now write x = I + w, f*x^j = -ln(1-I-w)*(I+w)^(j-1).

ln(1-I-w) = ln(1-I)+ ln(1 - w/(1-I)) = ln(1-I) - sum((w/(1-I))^k/k, k=1..infinity)

(I+w)^(j-1) = sum(binomial(j-1,k)*I^(j-1-k)*w^k, k=0..j-1)

...
 

> sys1:={diff(x1(t),t)=(3+5*I)*x1(t)-sqrt(2)*(cos(3)
  +I*sin(3))*y1(t)-2*I*z1(t),diff(y1(t),t)=conjugate(x1(t)),
  diff(z1(t),t)=0.25-0.3*z1(t)-2*I*(x1(t)-y1(t))};

I don't think dsolve likes that "conjugate".  Consider the real and imaginary parts.

> eval(sys1,{x1(t)=xr(t)+I*xi(t),y1(t)=yr(t)+I*yi(t),
     z1(t)=zr(t)+I*zi(t)});

> newsys:= map(e -> (map(evalc@Re,e), map(evalc@Im,e)), %);

> evalf(dsolve(evalf(newsys)));

Now for the original variables:

> [ x1(t)=eval(xr(t)+I*xi(t),%), 
    y1(t)=eval(yr(t)+I*yi(t),%), 
    z1(t)=eval(zr(t)+I*zi(t),%)];

Usually in complex analysis one would calculate the integral of a function over some particular contour or path.  You haven't said what the path is, and that's very important information in general.  In this case, why does the answer actually depend only on the endpoints of the path?  (That should be a big hint...)  And if the path is closed, it doesn't even depend on that...

For the histogram:

> Histogram(a,frequencyscale=absolute, discrete=true);
   H:= %:

For a smooth curve:

> plot(CurveFitting[Spline]([seq([i,numboccur(i,a)], i=0..35)], t),
      t=0..35);
   P:= %:

To display them together:

> plots[display](H,P);view=[-0.5 .. 35.5,default],axes=frame);

If you didn't have that "assume=integer", i.e. it really was a linear programming problem, I would suggest solving the dual problem.  A linear programming problem is unbounded if and only if it is feasible and its dual is infeasible.  The "assume=integer" makes it a bit trickier, but I think you can do this: a feasible
integer programming problem with rational coefficients is unbounded if and only if it is feasible and the corresponding homogeneous linear programming problem (i.e. with constant terms replaced by 0, and without the integer restriction) has feasible solutions with positive objective value.

For example, consider the problem:

maximize x1 
subject to
  -x1+2*x2 <= -1
  -3*x1+2*x2-3*x3 <= -8
  -4*x1-x2+x3 <= -13
  -3*x1+x2-2*x3 <= -9
   x1-4*x2+x3 <= -4 
  -x1-x2+4*x3 <= -1
   x1, x2, x3 >= 0

which I would tackle as follows:

> obj:= x1;
   cons:= {-x1+2*x2 <= -1, -3*x1+2*x2-3*x3 <= -8, 
     -4*x1-x2+x3 <= -13, -3*x1+x2-2*x3 <= -9, 
      x1-4*x2+x3 <= -4, -x1-x2+4*x3 <= -1};
   LPSolve(obj, cons, maximize=true, 
      assume={integer, nonnegative});

Warning, problem appears to be unbounded

           [6, [x1 = 6, x2 = 3, x3 = 0]]         

I don't believe the solution returned is "meaningless": in this case it is not a feasible solution, but it is a solution of the homogeneous problem with positive objective value.  But this does not seem to be documented, so maybe you don't want to trust it.

First look for a feasible solution, by making the objective 0.

> LPSolve(0, cons, assume=[integer, nonnegative]);

         [0, [x1 = 7, x2 = 3, x3 = 0]]

Now solve the homogeneous linear programming problem for objective >= 1:

> LPSolve(0, map(t -> (lhs(t) <= 0), cons) union {obj >= 1},
     assume=nonnegative);

                 [0., [x1 = 1.00000000000000022, x2 = .249999999741850580, x3 = 0.]]

Presumably that's an approximation to the actual solution x1 = 1, x2 = 1/4, x3 = 0.  The conclusion is that for any nonnegative integer n, we have a feasible integer solution [x1 = 7 + 4*n, x2 = 3 + n, x3 = 0]. 

This differential equation is unlikely to have any closed-form solutions.  DEtools[symgen] finds no symmetries.

Pts:= Array(0 ..100, 1..2, datatype= float[8]):
for i from 1 to 100 do 
  Pts[i]:= <i, Pts[i-1,2] + rand(-10..10)()>  
end do:
plots[pointplot](Pts,style=line);

> solve(sin(x)=0,x,AllSolutions);
                               Pi _Z1~

(where _Z1~ means an arbitrary integer).  Does the TI calculator really think it has to be a multiple of 3?

My guess is that py07ajpm wants to compare the solution of a system for one value of a parameter to the solution for another value.  The way to do this is to use two copies of the system (with different names for the variables) and use odeplot to plot the difference of a variable of one version with the variable in the other version.  For example:

sys:= {D(x)(t)=y(t)+lambda*x(t)*y(t),D(y)(t)=-x(t),x(0)=1,y(0)=0};
sys0:= subs(lambda=1.0 ,sys);
sys1:= subs(lambda=1.01, x=x1, y=y1, sys);
Sol:= dsolve(sys0 union sys1, numeric);
plots[odeplot](Sol, [t,x(t)-x1(t)], t = 0 .. 10);

It must be 1 + x + 2*x^2 + O(x^3), so the general solution (analytic in domain D) is F(x) = 1 + x + 2*x^2 + x^3*G(x)
where G(x) is an arbitrary function analytic in D.

(a) I don't think there is a difference.  As far as I can tell they look the same, and Maple treats them the same.

(b) I'm not sure what you mean.  Double quotes " are used as delimiters for strings.  While in TeX you can have left quotes and right quotes, in Maple there is only one ".

The convolution of f with itself is as follows:

> f:= x -> exp(3*x)/(1+exp(x)^5);
  ff:= int(f(x)*f(y-x),x=-infinity..infinity) assuming y::real;

    ff := y*exp(3*y)/(exp(5*y)-1)

This runs into numerical difficulties for y near 0 because of "catastrophic cancellation": when y is near 0, exp(5*y) is near 1, so the relative accuracy
of exp(5*y) - 1 will be poor.  For such values of y, it would be better to use e.g.
the following Maclaurin polynomial, which will be accurate to within 10^(-10) for
-0.03 <= y <= 0.03:

> ff7:= convert(series(ff, y, 7), polynom);

   ff7 := 1/5+1/10*y-11/60*y^2-1/10*y^3+91/720*y^4
      +43/600*y^5

Or you could use the following (produced with the help of minimax in the numapprox package) which is accurate to within 1.7*10^(-14) on the same interval:
 

  ffm := .199999999999999+.999999999962622e-1*y
      -.183333333318560*y^2-.999999667658748e-1*y^3
      +.126388820584991*y^4+.715927882186684e-1*y^5
      -.815549317855673e-1*y^6

 So take your convolution as

   piecewise(y < -0.03, ff, y < 0.03, ffm, ff)

and you shouldn't run into problems.

> DEtools[convertsys]({diff(x(t),t,t) + 3*diff(x(t),t)^2 
         + 6*x(t) = 1}, [],[x(t)],t,x,v);

                                  2
  [[v[1] = x[2], v[2] = 1 - 3 x[2]  - 6 x[1]],

                             d
        [x[1] = x(t), x[2] = -- x(t)], undefined, []]
                             dt

> subs(%[1], [Diff(x[1],t)=v[1],Diff(x[2],t)=v[2]]);

            d               d                   2
           [-- x[1] = x[2], -- x[2] = 1 - 3 x[2]  - 6 x[1]]
            dt              dt

You shouldn't use a formal parameter (which is what x and y are in this procedure) as a local variable. 

What happens when a procedure is called is that the formal parameters in the body of the procedure are replaced by the actual parameters.  In this case x is replaced by 6 and y is replaced by 10.  That results in a statement

6 := 10

which produces an error message.

If your actual parameter was something that could be assigned a value (e.g. an unassigned name such as X), there wouldn't be an error, and X would be assigned the value 10, but the print(x) statement would produce X rather than 10.  This is because of Maple's evaluation rules: normally, after the substitution occurs, the formal parameters are not evaluated again.  You could use eval to force evaluation.  But the simplest rule to avoid any trouble is: after assigning a value to a formal parameter, don't try to access the value of that parameter.   Instead, use a local variable to store any value that you will need to access.