You may be interested in a Taylor-series approach that works for nonlinear delay-differential equations. 

See this article

The equation of a cylinder of radius r centred on the line given parametrically by P + t U where U is a unit vector is ||R - P||^2 - ((R-P).U)^2 = r^2.

f:= (x,y,z) -> (x-p1)^2 + (y-p2)^2 + (z-p3)^2 -  
   ((x-p1)*u1+(y-p2)*u2+(z-p3)*u3)^2-r^2;

residuals:= [seq(f(data[i,1],data[i,2],data[i,3]),i=1..20)];

Optimization[LSSolve](residuals,{u1^2+u2^2+u3^2=1},r=0..20);

[6808.56703138011744, [r = 6.60008795660578507, p1 = 16.5125569207732036, p2 = 2.60785435102466145, p3 = -23.4751689199507184, u1 = .979678485137572274, u2 = -.182604346141747226, u3 = .829802297530015666e-1]]

But this is only a local optimum.  For a global optimum, you can use the Global Optimization Toolbox.

GlobalOptimization[GlobalSolve](add(residuals[i]^2,i=1..20),
  {u1^2+u2^2+u3^2=1}, p1=-20..20, p2=-20..20, p3=-20..20,
   u1=-1..1,u2=-1..1,u3=-1..1,r=0..20);

[2790.34734321769520, [p1 = -3.83047439918606525, p2 = -0.679459181163955538,   p3 = -20., r = 14.0306328235989746, u1 = -0.626422891743573240,   u2 = 0.0313602023358039925, u3 = 0.778852291456809452]]
 

To visualize it:

sol:= %[2];

P := subs(sol, <p1,p2,p3>);

U:= subs(sol, <u1,u2,u3>);

with(LinearAlgebra):

get an orthonormal basis U, V, W

V1 := CrossProduct(U,<1,0,0>);
V := V1/sqrt(V1[1]^2 + V1[2]^2 + V1[3]^2);
W:= CrossProduct(U,V);

parametric expression for the cylinder

cyl := P + s*U + subs(sol, r) * (cos(t)*V + sin(t)*W);

minimum and maximum values for s in the data

[min,max](seq((data[i,1..3] - P^%T) . U, i = 1..20)); 

Cyl := plot3d(cyl, s= floor(%[1]) .. ceil(%[2]),t=0..2*Pi):
with(plots): 
PP := pointplot3d(data, symbol=box, colour=black):
display([Cyl, PP], scaling=constrained);

I too would favour the solution using seq, but you could also do something like this:

> plots[listplot](map(sin, [$1..12]));

The exponential function is exp.  You don't need an extra ^.  So it should be

> J:=Int(exp(-a*x^2)*sech(x/2)^2, x=-infinity..infinity);

> value(J) assuming a > 0;

Maple doesn't seem to know this integral.  Do you have reason to think it has a closed-form expression?

This is a linear homogeneous system, so the solutions will form a linear space.  We can look for fundamental solutions, forming a basis of this linear space, in the form x(t)=A*exp(r*t), w(t)=B*exp(r*t), with A and B not both 0.  

dde:= {w(t)-(x(t)+D(x)(t)+w(t-a)),
    (D@@2)(x)(t)-D(x)(t)-w(t)};
eval(dde,{x=(t -> A*exp(-r*t)),w=(t -> B*exp(-r*t))});
expand(map(`*`,%,exp(r*t)));
 

{A*r^2+A*r-B, B-A+A*r-B*exp(r*a)}

From the first equation, A must be nonzero (and we may as well take it to be 1),
and then B = r^2+r.

eval(%, {A=1, B=r^2+r});
Q:= op(% minus {0});

Q := r^2+2*r-1-(r^2+r)*exp(r*a)

Now as a function of r for fixed nonzero a, this will have infinitely many roots.
For example, take a=1:

RootFinding[Analytic](eval(Q,a=1),r=-100 .. 100+100*I);

produces 17 roots, of which one is real; note also that the complex conjugate of a root is also a root. 

As t -> infinity, the dominant terms in a typical solution will be the ones for the roots with the greatest real part: these are .40661050935585+.51398193551780*I
and its complex conjugate.

Try collect(%, wsq).

I don't know about kakuro, but you might try the approach Joe Riel used for sudoku here

To handle it with Optimization, the problem would have to be expressed in terms of real variables, so you would consider real and imaginary parts of your complex variables. 

Find your maple11.ini file (which may be in the Users subdirectory of your Maple 11 directory) and change the PlotBGColor line from

PlotBGColor=default

to something like

PlotBGColor=0 111 222

(each of the three components should be an integer from 0 to 255).

This will affect plots to the "window" device in the Classic GUI; it does not affect inline plots, or plots in Standard GUI.  I haven't checked whether it affects other plot devices such as gif.

In other words, you want G to be an antiderivative of g with G(x) -> 1 as x -> infinity.

Let Q be the result of your simplify command.

1) Substitute D(G) for  g and INF for infinity (this is kind of cheating...).

  Q1 := eval(Q, {g = D(G), infinity=INF});

2) Tell int not to worry about discontinuities.

  Q2 := subsindets(Q1, specfunc(anything,int), 
      y -> int(op(y),continuous));

3) Tell Maple about the limit.

  Q3 := eval(Q2, G(INF) = 1);

Q3 := pi1-pi2*(-G(-(-q*w-q*pi2+pi2*C)/w)+1)

The animate command gives floating-point values to the animation variable.  If that's to be used as an index, this  would be bad.  However, you can use round to make an integer out of it (and use the frames option to ensure that the animation variable is incremented by 1 each time).  Thus:

pleaseWork := proc(x) pointplot3d([[p[round(x), 1],
  p[round(x), 2], x]], color = blue, symbol = box, 
  symbolsize = 40) end proc;

animate(pleaseWork, [x], x=1..20, frames=20);

Actually, the results of solve here are all real numbers, but their expressions in terms of radicals use complex numbers.  This fact, discovered in the 16th century by Italian mathematicians (notably Cardano, Tartaglia and Bombelli), was actually the main reason for the development of complex numbers.

If you want real expressions for your real roots, you can get them in terms of trig functions using evalc:

> evalc([solve(x^3-6*x^2-7*x+58=0)]);

[2/3*57^(1/2)*sin(1/3*arctan(1/126*4701^(1/2))+1/6*Pi)+2, -1/3*57^(1/2)*sin(1/3*arctan(1/126*4701^(1/2))+1/6*Pi)+2-1/3*3^(1/2)*57^(1/2)*sin(-1/3*arctan(1/126*4701^(1/2))+1/3*Pi), -1/3*57^(1/2)*sin(1/3*arctan(1/126*4701^(1/2))+1/6*Pi)+2+1/3*3^(1/2)*57^(1/2)*sin(-1/3*arctan(1/126*4701^(1/2))+1/3*Pi)]

This is a great example of the problem with implied multiplication.  Your expression copies into a 1D input region as

(r(r+(1/2)*s))(sin(t)/(r+(1/2)*s)+(diff(y(x), x))/(r^2+(1/2)*rs))

which shows three separate errors that are endemic to 2D input:

1) your first r is used as a function: r(...) rather than r * (...)

2) Similarly, (r(...))(sin...) is a function call rather than multiplication

3) rs is a single name rather than r*s

And none of the three errors is apparent until you look at the expression in 1D.

The cure in all three situations is to insert a space (or *): between r and (, between ) and (, and between r and s.

Implied multiplication doesn't really mean that you can omit the *, it just means that you can replace it with a space.  Well, in some situations you can get away with not using a space, but that just makes life more complicated.

Surely that's not what you meant. 

DiagonalMatrix(E) would not work, because DiagonalMatrix expects a list or Vector rather than a float.  Also MatrixAdd expects two Matrices, and (if the DiagonalMatrix command and the + worked) it's only given one.

I don't know what's causing the particular problem, except that the fact that the switch happens at Digits = 15 or so is probably significant: it means LinearAlgebra operations will be using software floats rather than hardware floats.  Perhaps there's some remember table or cache table that's filling up and not being cleared at garbage-collection time?

I don't find it slow at all.  For example:

Q:= z -> hypergeom([1],[1.5],-z);

L:= [seq(RandomTools[Generate](float(range=0.5e-6 .. 5e-6)), 
    count=1..10000)]:

ti:= time(): for t in L do Q(t) end do:  time()-ti;

         0.875

So on my computer it took 0.875 seconds of CPU time to do 10000 of these hypergeom evaluations.  I wouldn't consider that slow at all, as Maple functions go.  Compare that to the equivalent "closed form" version:

Q2:= z -> evalf(-1/2*Pi^(1/2)*(-1+erfc((-z)^(1/2)))/exp(z)
     /(-z)^(1/2));

ti:= time(): for t in L do Q2(t) end do:  time()-ti;

      12.453

I wonder if there's really something else going on.  Perhaps if you posted your actual code we could have something to suggest.