Mathematicaの強力な記号操作能力を使って多体問題を解くプログラムを作ります。
#粒子のリストを作成する
粒子は質量、電荷、位置、運動量の変数名、初期位置の変数名、初期位置、初期運動量の情報を持つとする。

##putindex
添字つきの変数を作成する関数putindexを作ります。

putindex[x_,i_] :=ToExpression[ToString[x] <> ToString[i]];

実行結果

##makeparticles
粒子のリストを作成する関数makeparticlesを作成する。Mapを使ってputindexをそれぞれの項目に作用させ、添字つきの変数名を作成している。Map[f,list,{-1}]はlistの最下層の要素にfを作用させるという意味です。

makeparticles[n_] :=
 Table[
  Map[
   putindex[#, i]&,
   {
    m, (*質量*)
    q, (*電荷*)
    {x, y}, (*位置*)
    {px, py}, (*運動量*)
    {x0, y0}, (*初期運動量*)
    {px0, py0} (*初期運動量*)
   }, {-1}
  ],{i, n}
];

実行結果

##makeparams
粒子の質量や初期位置の値を入れるルールのリストを作成する関数makeparamsを作成する。MapAt[f,list,{All,1}]はリストの第二レベルの1番目にfを作用させる。
Flattenは入れ子のリストを開封します。例:Flatten@{a,{b,c}}->{a,b,c} f@aはf[a]と同じ意味です。

makeparams[n_] := Flatten@Table[
 MapAt[
  putindex[#, i] &,
   {
    m -> RandomReal[{0.1, 1}],
    q -> RandomReal[{-0.1, 0.1}],
    x0 -> RandomReal[{-1, 1}],
    y0 -> RandomReal[{-1, 1}],
    px0 -> RandomReal[{-0.1, 0.1}],
    py0 -> RandomReal[{-0.1, 0.1}]
   }, {All, 1}
  ],
  {i, n}     
];

実行結果

->はRuleを表す記号で、Replace:/.を使ってこのルールを好きな式に作用させれます。例：f[a]/.a->b->f[b]。makeparamsで作ったルールをmakeparticlesで作ったリストに作用させるとさっき作ったリストのm1,q1などに数値が入ったリストができます。
最初から数値を入れていないのは、後でこの変数を組み合わせた式を作ったとき、記号式の状態で見れるようにするためです。こうするとデバッグなどで便利です。

実行結果

#ハミルトニアンを記述する。
ハミルトニアンとは系の全エネルギーを粒子の位置と運動量で表したものです。ハミルトニアンが記述できれば、正準方程式として機械的に解くべき微分方程式が得られます。
準備としてダサいですが、makeparticlesで作ったリストの何番目に何の項目が入ってるかを添字に直接名前を着けることで分かりやすくします。

mass = 1; charge = 2; position = 3; momentum = 4; position0 = 5;
momentum0 = 6;

##運動エネルギー関数T

T[particles_] := 
Sum[Total@
 ((particles[[i, momentum]]^2)/(2 particles[[i, mass]])),
 {i, Length@particles}
]

実行結果

##相互ポテンシャル関数V
ここでは静電力による相互作用を考えます。（プラズマを考える）

V[particles_,i_, j_] := 
(particles[[i, charge]]*particles[[j, charge]])
/Sqrt[Total@((particles[[i, position]] - particles[[j, position]])^2) + 0.01];

実行結果

粒子どうしが完全に重なったときにエネルギーが無限大になる問題を防ぐため、分数の分母のルートの中に0.01を足しています。こうすると穏やかなポテンシャルになります。

##ハミルトニアンH

H[particles_] :=
 T[particles] +
  Sum[If[i > j, V[particles,i,j], 0], 
  {i, Length@particles}, {j,Length@particles}
]

実行結果

#正準方程式を記述する。
まず位置と運動量の変数のリストを作ります。

ハミルトニアンもこれに合わせます。

さて、正準方程式はこういう形をしています。

\begin{align*}
\dot{q}_i &= \frac{\partial H}{\partial p_i} \\
\dot{p}_i &= -\frac{\partial H}{\partial q_i}
\end{align*}

まずこの式の左辺は、単純にpst全体を微分すれば作れます。右辺に関しては、まずハミルトニアンをpstで微分したリストを作り、これを半分に分割し、前後を入れ替え、位置による微分には全てマイナス1を掛けることで作れます。
これはMathematicaの強力なリスト操作関数群を使ってこう書けます。

Thread[
 D[pst, t] == 
  Flatten[Reverse@Partition[D[Ht, {pst}], Length@ps/2]*{1, -1}]
 ]

この部分を理解したい人は偏微分を通常表記にしてくれる関数(下記リンクのpdconv)を使ってこの式の各部分を調べてみてください。
http://blog.wolfram.com/2011/12/15/mathematica-qa-series-converting-to-conventional-mathematical-typesetting/

pdconv[f_] :=
 TraditionalForm[
  f /. Derivative[inds__][g_][vars__] :>
    Apply[Defer[D[g[vars], ##]] &, 
     Transpose[{{vars}, {inds}}] /.
      {{var_, 0} :> 
        Sequence[], {var_, 1} :> {var}}]
  ]

#`NDSolve`で微分方程式を解く。正準方程式を立式して微分方程式を解いてくれる関数`hamiltonSolve`を作ります。 Methodに関してですが、"SymplecticPartitionedRungeKutta"というのが正準方程式に特化したソルバらしいのですが、Mathematica 8ではバグ（下記URL）があるらしく使えませんでした。https://mathematica.stackexchange.com/questions/90066/symplecticpartitionedrungekutta-shows-strange-error?noredirect=1&lq=1

hamiltonSolve[H_, particles_, tcond_] :=
  
  Module[{ps, pst, ips, t, te, Ht}, {t, te} = tcond;
   ps = Flatten[{particles[[All, position]], 
      particles[[All, momentum]]}];
   ips = Flatten[{particles[[All, position0]], 
      particles[[All, momentum0]]}];
   pst = Through[ps[t]];
   Ht = H /. Thread[ps -> pst];
   NDSolve[
    Join[ 
     Thread[
      D[pst, t] == 
       Flatten[Reverse@Partition[D[Ht, {pst}], Length@ps/2]*{1, -1}]
      ],
     Thread[(pst /. t -> 0) == ips]
     ], ps, {t, 0, te},
    (*Method->{
    "SymplecticPartitionedRungeKutta","DifferenceOrder"->2,
    "PositionVariables"\[Rule]Drop[pst,Length@pst/2]}*)
    (*Method->{"ImplicitRungeKutta",
    "DifferenceOrder"->4},*)
    (*Method->{"ExplicitRungeKutta","DifferenceOrder"->4},*)
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 4},
    StartingStepSize -> 0.01,
    MaxStepSize -> 0.1,
    MaxSteps -> Infinity
    ]
   ];

実行結果

#可視化
電荷を色で、質量を円の大きさで表しました。

DrawParticles[particles_, sol_, t_] :=
 Graphics[
  {
   {
       Blend[{Red, Blue}, Rescale[#[[charge]], {-0.1, 0.1}]],
       Circle[
        {#[[position, 1]][t], #[[position, 2]][t]}, #[[mass]]
        ]
       } & /@ particles /. sol
   , Circle[{0, 0}, 10]
   }
  , PlotRange -> {{-15, 15}, {-15, 15}}
  ]

Manipulate[
 DrawParticles[particles /. params, sol, t],
 {t, 0, 100}
 ]

#結果
このままだと粒子が遠くにいってしまうので、原点から半径10の位置に急激に増加するポテンシャルV2をつけて壁にします。壁はTanh関数を使って作りました。

V2[particles_] := 
  Sum[Tanh[(Total@((particles[[i, position]])^2) - 100)/0.1], {i, 
    Length@particles}];

実行結果（10体問題）

青はマイナスの電荷、赤はプラスの電荷を表します。赤と青の丸がちゃんと引き合っているのが分かります。

読んでくれた方ありがとうございました。
参考ページ
https://mathematica.stackexchange.com/questions/102890/solve-motion-from-hamiltons-equations
http://eman-physics.net/analytic/hamilton.html
#コード全体

putindex[x_, i_] := ToExpression[ToString[x] <> ToString[i]]
makeparticles[n_] :=
  Table[
   Map[
    putindex[#, i] &,
    {
     m,
     q,
     {x, y},
     {px, py},
     {x0, y0},
     {px0, py0}
     }, {-1}
    ],
   {i, n}
   ];
makeparams[n_] := Flatten@Table[
    MapAt[
     putindex[#, i] &,
     {
      m -> RandomReal[{0.1, 1}],
      q -> RandomReal[{-0.1, 0.1}],
      x0 -> RandomReal[{-1, 1}],
      y0 -> RandomReal[{-1, 1}],
      px0 -> RandomReal[{-0.1, 0.1}],
      py0 -> RandomReal[{-0.1, 0.1}]
      }, {All, 1}
     ],
    {i, n}
   ];
mass = 1; charge = 2; position = 3; momentum = 4; position0 = 5; momentum0 = 6;
V[particles_, i_, j_] := (
  particles[[i, charge]]*particles[[j, charge]])/Sqrt[
  Total@((particles[[i, position]] - particles[[j, position]])^2) + 
   0.01];
V2[particles_] := 
  Sum[Tanh[(Total@((particles[[i, position]])^2) - 100)/0.1], {i, 
    Length@particles}];

T[particles_] := 
 Sum[Total@((particles[[i, momentum]]^2)/(
   2 particles[[i, mass]])), {i, Length@particles}]
H[particles_] :=
  T[particles] +
   Sum[If[i > j, V[particles, i, j], 0], {i, Length@particles}, {j, 
     Length@particles}] + V2[particles];
hamiltonSolve[H_, particles_, tcond_] :=
  
  Module[{ps, pst, ips, t, te, Ht}, {t, te} = tcond;
   ps = Flatten[{particles[[All, position]], 
      particles[[All, momentum]]}];
   ips = Flatten[{particles[[All, position0]], 
      particles[[All, momentum0]]}];
   pst = Through[ps[t]];
   Ht = H /. Thread[ps -> pst];
   NDSolve[
    Join[ 
     Thread[
      D[pst, t] == 
       Flatten[Reverse@Partition[D[Ht, {pst}], Length@ps/2]*{1, -1}]
      ],
     Thread[(pst /. t -> 0) == ips]
     ], ps, {t, 0, te},
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 4},
    StartingStepSize -> 0.01,
    MaxStepSize -> 0.1,
    MaxSteps -> Infinity
    ]
   ];
DrawParticles[particles_, sol_, t_] :=
  Graphics[
   {
    {
        Blend[{Red, Blue}, Rescale[#[[charge]], {-0.1, 0.1}]],
        Circle[
         {#[[position, 1]][t], #[[position, 2]][t]}, #[[mass]]
         ]
        } & /@ particles /. sol
    , Circle[{0, 0}, 10]
    }
   , PlotRange -> {{-15, 15}, {-15, 15}}
   ];

particles = makeparticles[10];
params = makeparams[10];
sol = hamiltonSolve[H[particles] /. params, particles /. params, {t, 100}]
Manipulate[
 DrawParticles[particles /. params, sol, t],
 {t, 0, 100}
]