In principle the Statistics package should be able to do this sort of thing. However, with symbolic parameters s11,...,s22 I suspect this may be very complicated to calculate. There are too many special cases. Also, have you considered that your Z is not always going to be real (e.g. if y1 < 0)?

Left quotes `` are needed around `not` because not is a reserved word. The left quotes make it into a name that can be manipulated like other functions. @ is the composition operator, i.e. f @ g is the function such that (f @ g)(x) = f(g(x)). So `not` @ has is a function: (`not` @ has)(f,x) returns true if has(f,x) returns false and vice versa. A more pedestrian way of doing the same thing would be > SPLIT := (f,x) -> [selectremove(g -> not has(g,x), f)];

Note that the imaginary parts are very small (less than 10^(-6)). The equation leads to a cubic equation in sqrt(-18+6*F), and solving a cubic in terms of radicals uses complex numbers. The imaginary parts should really cancel out exactly, but in calculations with floats there are roundoff errors.

Once r has been assigned a value, it does indeed keep that value unless you explicitly remove it (e.g. by r := 'r', or by restart). I'd have to look at your worksheet to see where that is happening. But certainly it doesn't happen spontaneously.

In this type of problem, there's really no need for Maple, nor are there any special facilities in Maple for it. You have to formulate the problem in terms of probabilities and conditional probabilities. The main challenge is to be clear on exactly what are the rules governing the host's behaviour. You could use Maple as a calculator, putting the numbers into the formula and having Maple compute the answer, but in this case the calculation should be simple enough to do in your head.

My fault: I should have used

 print(display([WM], scaling=constrained));

(and then you can use the colon). I had done it first with a semicolon, and then tried to "clean it up" with colon, stupidly forgetting that this suppressed the printing of the display as well as the unwanted output of the other commands.

You can also use a mod b to get the remainder. Yes, iquo(a,b) is much more efficient than floor(a/b). Consider:

> tens:= [seq](10^m,m=1..10000): 
  sevens:= [seq](7^m,m=1..10000):
  ti:= time():
  for m from 1 to 10000 do floor(tens[m]/sevens[m]) od:
  time()-ti;

20.968

> restart;
  tens:= [seq](10^m,m=1..10000): 
  sevens:= [seq](7^m,m=1..10000):
  ti:= time():
  for m from 1 to 10000 do iquo(tens[m], sevens[m]) od:
  time()-ti;

1.016 On the other hand, there seems to be little difference in efficiency between mod and irem.

>  restart;
  tens:= [seq](10^m,m=1..10000): 
  sevens:= [seq](7^m,m=1..10000):
  ti:= time():
  for m from 1 to 10000 do tens[m] mod sevens[m] od:
  time()-ti;

0.999

>  restart;
  tens:= [seq](10^m,m=1..10000): 
  sevens:= [seq](7^m,m=1..10000):
  ti:= time():
  for m from 1 to 10000 do irem(tens[m],sevens[m]) od:
  time()-ti;

0.985

For each date you want the plot of the world map together with the circles for all dates up to and including that one? And may I assume that the dates in your file are in order (as they are in your example)? Then you can try this. Your Array is E.

> with(plottools): with(plots):
  Globe:= (lon,lat) ->(cos(lon)*cos(lat),  
      sin(lon)*cos(lat), sin(lat));
  deg:= Pi/180;
  circonsphere := ... as defined earlier ...
  WM:= WorldMap_Complete:
  N:= rhs(rtable_dims(E)[1]):
  i:= 1: 
  while i <= N do
    date:= E[i,1];
    for j from i+1 to N do
      if E[j,1] <> E[i,1] then break end if
    end do;
    WM:= WM, seq(
       circonsphere([E[k,2]*deg,E[k,3]*deg], 0.05, 
        color=red, thickness=1), k=i..j-1);
    plotsetup(gif, plotoutput=cat("WorldMap_Test2A_",
      StringTools[SubstituteAll](date,"/","-"),".gif"),    
      plotoptions="height=300,width=600");
    display([WM], scaling=constrained);
    i:= j;
  end do:
  plotsetup(default);

You might try the book "Introduction to Vectors and Tensors, Vol. 1, Linear and Multilinear Algebra" by Ray M. Bowen and C.-C. Wang, which can be downloaded for free at <http://www1.mengr.tamu.edu/rbowen/>

Instead of using the Tutor interface, you might try Hint. For example:

> with(Student[Calculus1]):
  J1 := Int(tan(x)*sin(x),x);

J1 := Int(tan(x)*sin(x), x)

> Hint(J1);

[change, u = cos(x), u]

> J2 := Rule[%](J1);

J2 := Int(tan(x)*sin(x), x) = Int(-sqrt(1-u^2)/u, u)

This system of ODE's should be OK except at (x,y) = (0,0). The first equation (v[0]^2*y(t)^2/(-diff(y(t),t)+w)^2)^(1/2) = v[0]*y(t)/(-diff(y(t),t)+w) will be true as long as the right side is positive, and then the other equations determine x(t) in terms of y(t). So it seems to me that there's a bug in dsolve here. Note that if we write x(t) = r(t)*cos(theta(t)) and y(t) = r(t)*sin(theta(t)), we get what seems to be a correct result.

> polareqs:= simplify(eval({sys_ode}, 
   {x(t) = r(t)*cos(theta(t)),
    y(t) = r(t)*sin(theta(t))}),symbolic);

polareqs := {diff(r(t),t)*cos(theta(t))-r(t)*sin(theta(t))*diff(theta(t),t) = -v[0]*cos(theta(t)), diff(r(t),t)*sin(theta(t))+r(t)*cos(theta(t))*diff(theta(t),t) = -v[0]*sin(theta(t))+w}

> S := simplify(dsolve(polareqs union {r(0)=a,  
   theta(0)=0}));

S := {r(t) = a/cos(RootOf(-Int(1/cos(_a)^2*((1+sin(_a))/cos(_a))^(-1/w*v[0]),_a = 0 .. _Z)*a+w*t))*(1/cos(RootOf(-Int(1/cos(_a)^2*((1+sin(_a))/cos(_a))^(-1/w*v[0]),_a = 0 .. _Z)*a+w*t))*(1+sin(RootOf(-Int(1/cos(_a)^2*((1+sin(_a))/cos(_a))^(-1/w*v[0]),_a = 0 .. _Z)*a+w*t))))^(-1/w*v[0]), theta(t) = RootOf(-Int(1/cos(_a)^2*((1+sin(_a))/cos(_a))^(-1/w*v[0]),_a = 0 .. _Z)*a+w*t)}

> odetest(S, polareqs);

{0}

Drag-and-drop doesn't leave a record of what second function was plotted (i.e. the worksheet has only the command that produced the first plot), and re-executing the worksheet will only give you the first plot. So I much prefer to do it the old-fashioned way:

> plot([sin(x), cos(x)], x=0..2*Pi);

You also have easy control over what colours are used:

> plot([sin(x), cos(x)], x=0..2*Pi, colour=[red,blue]);

What could be simpler?

When given a system, dsolve returns a set. A set has no fixed order, and you shouldn't write code that depends on it being in a certain order. What you can do instead is construct a list, e.g. like this:

> Set:= dsolve(whatever);
  List:= [F(x) = subs(Set, F(x)), G(x) = subs(Set, G(x))];

As for the constants, I think the best way of fixing them is to specify initial conditions (using symbolic constants) as well as differential equations, e.g. F(0) = C1, D(F)(0) = C2, ....

I presume the entries of M are polynomials in the variables k1, k2, omega. Then N is a (probably rather complicated) polynomial in those variables, and probably can't be solved in closed form. If a plot is what you're after, I'd suggest either using implicitplot or fsolve. Something like this:

> H:= subs(k1=0, k2=k, N): 
  plots[implicitplot](H, k = 0 .. 0.5, omega = -100 .. 100);
# or whatever interval for omega you think appropriate

or

> H:= subs(k1=0, k2=k, N):
  kvals:= [seq](j*0.5/100, j=0..100):
  pts:= map(t -> seq([t,y],y=[fsolve(subs(k=t,H))]), kvals):
  plots[pointplot](pts, symbol=point);

This is indeed the equation of a cone, or more precisely a pair of cones touching at the vertex. But it will look more cone-like if you change the bounds so the cones are cut off on planes x=constant. Note that for x = 50, say, y will go from -200 to 200 and z from -100 to 100. It also helps if you use a fairly large grid setting and make sure the origin is a grid point. Try this:

> with(plots):
 implicitplot3d(16*x^2=y^2+4*z^2, x=-50..50, y=-200..200,
  z=-100..100, grid=[31,31,31], style=patchnogrid, axes=box);