isolve is rather limited in its capabilities. This is not entirely its fault, because there is no algorithm to find solutions to diophantine equations (look up Hilbert's 10th problem). But there are even some two-variable quadratic cases that isolve can't do (though they should be reducible to Pell's equation). For example:

> isolve(x^2-x*y-y^2=1);

For the Hessian matrix, use Hessian in the VectorCalculus package. Partial derivatives of any order can be found using diff.

Do you mean something like this?

> plots[pointplot]([seq(seq(seq([f[i,j,k],t[i,j,k]],
      i=-n..n),j=-n..n),k=-n..n)]);

It means that you consider the coefficients of the polynomial as members of the field Z[13] of integers mod 13. In more concrete terms, you are looking for polynomials whose product, after reducing the coefficients mod 13, is the given polynomial. As Scott said, you can use Factors(...) mod 13. You can also use Factor(...) mod 13, which will return a product of factors rather than a list with factors and multiplicities.

You could try something like this.

> f := proc(a,b)
   local ca, cb, aa, bb;
   if type(a,`+`) then return map(procname,a,b)
   elif type(b,`+`) then return map2(procname,a,b)
   end if;   
   if type(a,`*`) then
     ca,aa:= selectremove(type, a, numeric);
   else
     ca:= 1; aa:= a
   end if; 
   if type(b,`*`) then
     cb,bb:= selectremove(type, b, numeric);
   else
     cb:= 1; bb:= b
   end if;
   ca*cb*f2(aa*bb);
 end proc;
> f2:= proc(ab)
    option remember;
    int(int(ab,t=0..tf),tf=0..T);
  end proc;    

Depending on what type of coefficients you anticipate, you might replace numeric with freeof(t).

If the assumptions to be made will change from one command to the next, it may be better to use "assuming" rather than "assume". It may be useful to use a variable to keep track of the current assumptions. For example:

> assumptions:= A >= 0, B > 0;
  sqrt(B^2) assuming assumptions;

> for i from 4 to 1 by -1 do
   i
  end do;

LinearAlgebra does make extensive use of compiled external code, the NAG routines. In fact, this might be connected to the troubles you're having in profiling. See the help page "Efficient Numeric Linear Algebra Computation" (?EffNumLA) for pointers on how to make your code more efficient. In particular, it may make a significant difference if you use datatype=float[8].

I agree with Scott that it would be better to post a worksheet, but I'll take a stab at it anyway. It looks to me like the solution of your equation becomes undefined (or more precisely, Maple's numerical approximation of the solution becomes undefined) at about t = 0.2. Why that could be, I can't say without looking at your equation.

After using sin(e/lambda)^2 = 1/2 - 1/2*cos(c/lambda) and substituting sqrt(1-k[2]^2) = r, I find that the integral in question is a constant multiple of Int((1-alpha*r*cos(c/lambda)+beta*r^2)/(1+gamma*r^2-delta*r*cos(c/lambda)),lambda ) with alpha = 2.494476546, beta = 1.555603312, gamma = .4046124214, delta = 1.272183039. The change of variables c/lambda=t makes the definite integral into a constant multiple of Int((1-alpha*r*cos(t)+beta*r^2)/(1+gamma*r^2-delta*r*cos(t))/t^2,t = 38446.38339 .. 38751.51341) Note that the endpoints are large and relatively close together (though their difference is nearly 100 pi). So over this interval the factor 1/t^2 is nearly constant, and maximizing the integral should be mainly a matter of maximizing S = Int((1-alpha*r*cos(t)+beta*r^2)/(1+gamma*r^2-delta*r*cos(t)),t=0..2*Pi) Now Int(1/(1-k*cos(t)),t=0..2*Pi) = 2*Pi/sqrt(1-k^2) for abs(k) < 1 so I get S = -2*Pi*(-alpha*(1+2*gamma*r^2-delta^2*r^2+gamma^2*r^4)^(1/2)+alpha+alpha*gamma*r^2-delta-delta*beta*r^2)/delta/(1+2*gamma*r^2-delta^2*r^2+gamma^2*r^4)^(1/2) Differentiating that with respect to r, I get r=0 and (+/-)1/(alpha*gamma*delta-delta^2*beta+2*beta*gamma-2*gamma^2)*((alpha*gamma*delta-delta^2*beta+2*beta*gamma-2*gamma^2)*(-delta^2-2*beta+alpha*delta+2*gamma))^(1/2) which would be outside the interval [0,1]. So I suspect that the maximum of P[through] is either at k[2] = 0 or 1.

The integral int((a+b*cos(c/lambda))/(1+d*sin(e/lambda)^2), lambda) has almost certainly no closed-form solution in general. It looks to me in Xan's case c = 2*e. I think there may be relations between the other constants as well, e.g. 1.428285686 is presumably meant to be 2 * .7141428429. If you want to get closed-form solutions, there's a better chance if you leave the constants in symbolic form rather than using floating point. Not that it's very likely that there's a closed form solution. Still, if you could do reliable numerical evaluations of the integral you should be able to use numerical methods to find the maximum. But the integrand is oscillating pretty wildly in this interval, so numerical evaluation is going to be hard. A higher setting of Digits is almost certainly needed; in fact, it may already be too late as the result could be quite sensitive to the roundoff error that has already occurred.

The real basin is a cross-section of the complex basin. Complex dynamics has some very interesting theory behind it. The boundary of the complex basin of attraction is often a nice interesting-looking fractal; a subset of the real line is not so interesting to look at. Anyway, you might look at something like this. Given your polynomial, construct a function Newt that performs a step of Newton's method starting at x. Plot some iterate of Newt, say Newt@@10, over some interval (using an interval for the y axis as well to avoid very large values) and using discont=true. The regions where the graph is almost horizontal at a y value which is a root will be in the basin of attraction of that root. For example:

> f:= x -> x^3 + 2*x^2 - x -2;
  Newt:= unapply(x - f(x)/D(f)(x), x);
  plot(Newt@@10, -3 .. 2, -3 .. 2, discont=true);

First you have to define precisely what you mean by "correlated random numbers from an arbitrary distribution". What joint distribution do you have in mind?

Hints: The command is called sum, but you don't have the right syntax. Look at the help page. For the second one, I suspect you're missing some parentheses.

To get piecewise linear functions from those lists:

> f:= unapply(CurveFitting[Spline](
      [[1, 0], [3, 3], [5, 3], [7, 0]],x,degree=1),x);
  g:= unapply(CurveFitting[Spline](
      [[1, 0], [5, 4], [7, 4], [9, 0]],x,degree=1),x);

To find the intersection:

> solve(f(x) = g(x));

1, 4, 15