Two orbiting planets in perpendicular planes Inspired by this question. Can a 3 body problem, starting with two planets orbiting a larger one (so massive it may be taken to stand still) in perpendicular planes, be stable?
Is there known an analytical solution to this 3 body problem?
Or a qualitative description of the evolution?
Will the two bodies approach coplanar orbits?  

Bounty for this:  
Is Classical Mechanics (that is Newtons 1/r^2 law) sufficient for some large random collection of point-particles (with nonzero net angular momentum) in orbit around a larger one, sufficient to explain that they approach coplanar orbits?
Or is it needed to add collisions and energy loss etc..
 A: Poincaré 

The version finally printed contained
  many important ideas which lead to the
  theory of chaos. The problem as stated
  originally was finally solved by Karl
  F. Sundman for n = 3 in 1912 and was
  generalised to the case of n > 3
  bodies by Qiudong Wang in the 1990s.

Karl_F._Sundman

used analytic methods to prove the
  existence of a convergent infinite
  series solution to the three-body
  problem in 1906 and 1909.

Qiudong_Wang

Wang is best known for his paper The
  global solution of the n-body
  problem (*), in which he generalised
  Karl F. Sundman's results from 1912 to
  a system of more than three bodies.

(*)With Zero Angular Momentum, it seems.
There are a large colection of N-Body codes available from the net, and some of them work with GPUs (graphics hardware)
a SoftPedia list of opensource codes 
I've downloaded Gravit from the site of Gerald Kaszuba:  
I've choosed his work because it is loaded with options, even if it is NOT physically correct:
I've included the Velocity Verlet Integrator to solve this problem.  
// Velocity Verlet integrator algorithm
// r(i+1) = r(i) + v(i)*dt + a(i)*dt^2 / 2
// v(i+1/2) = v(i) + a(i) * dt/2
//     1st -- act. pos, 2nd calc accel = a(i+1), 3rd vel, 4th acc
// a(i+1= sum of accels on i+1  (accel)
// v(i+1) = v(i+1/2) + a(i+1)*dt/2  

I've also included code to work with my GPU using BROOK+ of AMD/ATI (320 parallel processors, 60000 objects)
In the site of CUDA/OpenCL of NVIDEA can be found NBODY code to work in OpenCL.
Sverre Aarseth codes are a strong reference (code in Git).
Theory in The Art of Computational Science and Maya Project
Grav-Sim
types of code, and links resources 
There are solutions stables with 3 body as we can see in this image. 
Your question with two bodys in perpendicular orbits, I think that in general the orbits will tend to become planar.  I didnt tested yet about the stability but I think that for aproximate radius they will be unstable.  
EDIT add:
'The three-body problem with close encounters is notoriously ill-conditioned because it admits chaotic solutions that manifest extreme sensitivity to initial conditions.'
My code is good for a large ensemble of bodies like a galaxy, a star cluster or a disk of matter.
The code for a few bodies must be of 'direct' type (and non GPU) and make use of 'adaptive time step'. In the short range distances the time step must be finer then in the long range and the use of a different reference frame (a local) could benefit the solution.  
Picard and Parker-Sochacki based Methods
If I were you I'd explore this very interesting paper An adaptive N-body algorithm of optimal order (2002 ) and this one and this one and THIS ONE with 47 lines of code inside
(I must say 'thank you' because you made me find this new, interesting way) 

Table I: Basic Source Code for Solving
  the N-body Problem
  The preceding 47 lines of code ...
  What is stunning is
  the simplicity of the solution.

EDIT add 2
I've done some simulations available here where a xls can be downloaded.
 
A: If you have two orbiting bodies, with angular momentum ${\vec L}_1$ and ${\vec L}_2$ the total angular momentum of the system is ${\vec J}~=~{\vec L}_1~+~{\vec L}_2$.  The angular momentum of the system has a pseudo-potential $C(r_1,r_2){\vec J}\cdot{\vec J}$, for $C(r_1,r_2)$ a function of the orbital radii and masses of two bodies.  The relevant portion is the angular momentum squared which is
$$
{\vec J}\cdot {\vec J}=~{\vec L}_1\cdot {\vec L}_1~+~{\vec L}_2\cdot {\vec L}_2~+~2{\vec L}_1\cdot {\vec L}_2
$$ 
This is related to the orbit-orbit interaction in QM, but we can physically see that the total angular momentum $\vec J$ is largest when the two ${\vec L}_i$ are parallel, and minimal when they are anti-parallel.  The potential, or pseudo-potential term is negative, and so the potential is minimal (most negative) when the orbits are aligned with each other.  This is the opposite of what happens in QM with electrons, due to the magnetic poles associated with orbital and spin directions.
This is somewhat descriptive, for a full presentation requires working with the 3-body problem in classical mechanics.  That would require a much more extensive coverage of the problem.  However, this roughly answers the basic question.  It is also one reason that solar systems (stellar systems) have planets which orbit in approximately the same plane and with angular momentum approximately aligned.
A: Strictly speaking, no, if you are asking whether the two smaller planets can stay in orbits in perpendicular planes.  Looking down on the system along the line where the two planes intersect, at some stage the two will have an angle at the the megaplanet less than 180 degrees or $\pi$ radians between them, and at that point the gravitational attraction between the small planets will start to move them out of their original planes.
This does not preclude a more complicated pattern where the orbits of the two smaller planets are not in planes but are something more complicated with substantially different average inclinations relative to each other.  I don't know, but perhaps it is possible for these to be perpedicular on average.  
A: There are free programs for playing with gravitational simulations  for example satellite orbits , so there is no problem of calculating numerically any kind of  orbits.
I do not see why two perpendicular orbit planes would not be stable in principle. It would depend on the parameters of course, masses and distances. If the satellite  bodies or the central one have oceans, tides will affect the orbits and long term stability will be lost. Otherwise, once the orbits are set up angular momentum conservation should keep them in the plane.
If you are really interested you could try downloading and playing around with some program to make sure.
Edit/correcting : in this wiki article the three body problem is examined and it is stated that the orbits are in general chaotic, in fact the three body gravitational problem led to the study of deterministic chaos. In the gif provided for a three  small mass particle problem I cannot see what goes on with angular momentum conservation. Some stability exists in special solutions.
So the problem is not as simple as observing the sun moon and earth, and depends on specific initial assumptions. 
A: In the limit of infinitesimal secondary masses the problem is clearly stable.
So consider the problem of how the secondary (planetary) masses interact with each other. For simplicity, let's make the situation symmetric about the plane midway (45 degrees) between the two planetary orbits. To minimize their interaction, let's make them 180 degrees out of phase so that one is at apogee while the other is at perigee.
If we treat the interaction between the planets as a perturbation, the system is stable before we include the perturbation (of course). Since the orbital periods are the same, the perturbation is stable. The question is "can we arrange for the perturbation to be zero?"
This is an interesting problem and one that it appears I will have to simulate in Java.
