Your final do-loop assigns columns, in-place, into what you intend to be Matrix P. But you haven't created the Matrix P, yet. Try adding,
P := Matrix(3,3):
before that final do-loop.
You might also think that P should be purely real. You could try simplify(map(fnormal,P)) at the end. But maybe you would rather compute the eigenvectors of a floating-point symmetric version of stif2, instead of doing floating-point evaluation of the eigenvectors of the exact stif2? Ie,,
Eigenvectors(evalf(Matrix(stif2,shape=symmetric)));
acer

You posted this in the category, "How do I ...? with Maple". So here's how to do it in Maple,
y := ln(sin(x)) - (1/2*sin(x)^2);
diff(y,x);
acer

Enter any of these at the maple prompt, and then press the key.
?scatterplot
?Statistics,scatterplot
?stats,scatterplot
acer

Sorry for not being clear enough

here. What I was trying to get across was that there may be neither a name nor some known, especially nice, closed form for the convergent result.
Just don't

name it after yourself. :)
acer

Maple 11.00, 64 Linux,
> ba,bu,st:=kernelopts(bytesalloc),kernelopts(bytesused),time():
> N:= 10^6:
> R := [ add( rand(0..1), k=1..12 )-6$N ]():
> kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu,time()-st;
376107760, 664232360, 39.998
New session (not restart),
> ba,bu,st:=kernelopts(bytesalloc),kernelopts(bytesused),time():
> N := 10^6:
> g := rand(0..1):
> R := Vector[row](N,datatype=integer[4]):
> for i from 1 to N do
> R[i] := add(g(),j=1..12) - 6;
> od:
> kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu,time()-st;
10352792, 217286880, 12.098
So, my posted method seems to use 1/36th the allocated memory, and take 1/3 of the time, for N=10^6.
There may be even leaner ways, perhaps using Statistics.
acer

Could it be that you don't really want to keep each of the one million sequences?
N := 10^6:
g := rand(0..1):
R := Vector[row](N,datatype=integer[4]):
for i from 1 to N do
R[i] := add(g(),j=1..12) - 6;
od:
After the above, the results are all stored in Vector R.
acer

There are various ways to accomplish this, and for your purposes you might be free to pick and choose.
I like the way that Joe suggested. But perhaps it's not to your taste, because you originally considered, map(x->1/x,test) rather than, say, map(`/`,test).
You might also wish to consider map(t->map(x->1/x,t),test) , if it's more to your style or liking.
acer

What value are you hoping to see instead, for the derivative at the point(s) where the functions are not continuous?
acer

Please realize, that you are not the first person to run amok of this issue.
By default, a subscripted name is (under the hood) a table reference. So your r_2 is actually r[2] underneath the pretty typesetting of 2D Math. But you generally cannot use the name r for one thing, and also use it as a table by referencing the indexed names like r[2]. It's a clash, trying to use the same base name in two ways.
Ok, so the name and its indexed form cannot generally be used at the same time. But you want to use the name, and its subscripted form, at the same time. So the question is, can subscripted names be made so that they are not indexed name forms underneath. The answer is yes, but it isn't as easy as it should be.
You want to use both the name r and the subscripted r names such as r_2, and not have them collide. One way to bring this about is to toggle the subscripted name as its own distinct name, rather than as an indexed form r[2] of the name r. You can use the select-and-right-click context menu on the subscripted typeset name, and turn it into its own uniquely identifying name (wholly distinct from, say, name r).
The context-menu action to do that is 2D Math -> Convert To -> Atomic Identifier.
That's not a very nice workaround, because it forces use of the mouse, makes you go about 3 levels deep in menus. And most unfortunate of all, it is necessitated that this action be repeated on all instances of r_2, say, everywhere in the document.
Another workaround might be to use the so-called subliteral entry in the Layout palette. It's hardly less work than the context-menu toggle.
This comes up so often here on mapleprimes. Clearly the default behaviour of Maple could be revisited. Maybe the default could be to have subscripts create unique names. Or some global control switch might be provided. Or some side pane might show which variables are in use in which way. And so on.
acer

p :=proc(N)
global X, Y;
local i,your_pick,prize,unchosen_doors,montys_door,
other_door,noswitchtotal,switchtotal,S,P,Two;
noswitchtotal,switchtotal := 0,0;
S := Statistics[Sample](X,N);
P := Statistics[Sample](X,N);
Two := Statistics[Sample](Y,N);
for i from 1 to N do
your_pick := trunc(P[i]);
prize := trunc(S[i]);
unpicked_doors := {1,2,3} minus {your_pick};
montys_choices := {1,2,3} minus {your_pick} minus {prize};
if nops(montys_choices) = 1 then
# Case where Monty has no choice (Player picked a loser)
montys_door := op(montys_choices);
else
# Case where Monty has a choice, and chooses at random.
# Sort his two choices, and pick 1st or 2nd at random.
montys_door := sort(convert(montys_choices,list))[trunc(Two[i])];
end if;
other_door := op(unpicked_doors minus {montys_door});
if prize = your_pick then noswitchtotal := noswitchtotal + 1; end if;
if prize = other_door then switchtotal := switchtotal + 1; end if;
end do;
noswitchtotal, switchtotal;
end proc:
X := Statistics[RandomVariable](DiscreteUniform(1,3)):
Y := Statistics[RandomVariable](DiscreteUniform(1,2)):
p(6666);
# Not switching wins 1/3 or the time, on average.
# Switching wins 2/3 of the time, on average.
acer

I should qualify my off-the-cuff answer. Programming simple things in a rich language like maple is often straightforward. But in this case, if one is not careful enough, then the program might simulate the "wrong" game, leading to the wrong conclusion about switching.
Someone might post such a program, but you'll have to give them a little while to write it.
Note that your original statement of the game is not precise, or usual, and may represent a "wrong" game here. It is usually stated more somethign like this: The player gets to pick any of the three doors. Monty then shows him behind one of the remaining doors, and that door will not have the prize. The player then has the choice to stick with the first choice, or to switch to the other unopened door.
I remember back in grad school, that one or two colleagues were obsessed with this game. A lot of the dispute centered on the exact description.
Not so long ago I saw it come up in practice, I believe. It was in one of Zia's newspaper columns on Bridge. South had to make some choices (finesses, say). The good play was to choose one path of three, then put the opponent on play in such a way that he was forced to disclose one of the others. The idea was then that the odds were better to then switch to the other plan.
acer

In Maple 11.00, this worked for me,
textplot([0.6,0.8,typeset(`#msubsup(mi("p"),mi("gt"),mo("*"))`)]);
Now, you might well ask how one produces such a beast as `#msubsup(mi("p"),mi("gt"),mo("*"))` in the first place.
I used the Layout palette, and was in 2D Maple notation mode (in a Worksheet, but it probably works in a Document too).
I selected the first entry in that palette, which looks like A raised to the b.
For A, I entered, p_gt , and then I hit tab to get to fill in the b. For b, I entered, \* . Then I hit return. That gave me the nice looking 2D input object.
But I found that I could not enter that, by simply copy 'n pasting the 2D object, in the typeset() call within textplot.
So instead I selected all of that nice looking 2D object, and used the 2-D Math -> Convert To -> Atomic Identifier choice in its context-menu.
Once having toggled it as an atomic name, I cut and pasted it into another blank line. I then used the context-menu on that, to convert it to 1-D maple notation input. That is to say, 2-D Math -> Convert To -> 1-D Math Input was the context-menu choice. The result was that funny piece of MathML-wannabe.
Hope that helps some.
acer

Here's how I see what you've tried.
You want to solve the position of the rocket, {X3,Y3,Z3}, out to 1 year's worth of 1-seconds (ie. problem is scaled so that independent variable t is in seconds).
The derivatives of the rocket's motion/position w.r.t will depend on the position of the two planets, {X,Y,Z} and {X1,Y1,Z1}.
Instead of building one large differential system to encapsulate all three objects (where only "influence" is of two upon the third) you instead solve each of the two planets' positions separately, using listprocedure output. Then you embed those resulting on-the-fly calls (to solve the planets' positions) into the differential equations for the rocket's motion/position. Doesn't that force maple to do more work than if the three systems were rolled up together?
You're using dsolve/numeric method=lsode throughout, a hard working (if somewhat antiquated) stiff solver. Is that the best idea here? (You may want routines for dX/dt, dX1/dt, etc, depending on how you go here.)
acer

Perhaps it was 16*x = y^2+4*z^2 that was meant instead?
acer

Vo/Vs isn't (technically) a subexpression of the original equation. But you can substitute for something like it, then isolate, and then substitute back again.
eval(normal(isolate(eval(Vs/R1=(Vs-Vo)/R2,Vo=Z*Vs),Z)),Z=Vo/Vs);
acer