下のソースプログラムをそのままテキストエディタ（メモ帳など）にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(この class file name は hyion なので、このテキストエディタを "hyion.java" とセーブしてコンパイルしてほしい。)
ここでは 1 MM = 1 × 10-14 meter という新しい単位を使用している。
このサンプルプログラムでは、実行すると、最初に H2+ の核間距離 ( MM ) を入力する。
その後、電子１の ｙ座標 "r" ( MM )、x 座標 "a" ( MM ) を 入力する。
これら入力値から、このプログラムは H2 の結合エネルギー (eV)、力 F1 と遠心力の比、n2 に作用する 電子と n1 からの力の比 を出力する。
F1 は 電子に作用する 原子核 1 (n1) 方向の力である。
電子に作用する力と その速度から １軌道に含まれるドブロイ波を計算する。
( ここでは ビリアル定理 E = -T = 1/2 V を使用している。)
さらに 原子核 2 の力による x 軸周囲の歳差運動を考慮して このドブロイ波を計算する。
total waves は これらの合計である。
また 原子核 2 からの x 方向の力の影響を考慮した ドブロイ波も計算する (= new total waves )
import java.util.Scanner;
class hyion {
public static void main(String[] args) {
Scanner stdIn=new Scanner(System.in);
System.out.println("nuc (MM) ? (= internuclear distance of H2+ ion");
double nnuc=stdIn.nextDouble();
System.out.println("r (MM)？ = electron y coordinate ");
double rr=stdIn.nextDouble();
System.out.println("a (MM) ? = electron x coordinate ");
double aa=stdIn.nextDouble();
double me=9.1093826e-31;
double pai=3.141592653589793; double epsi=8.85418781787346e-12;
double h=6.62606896e-34; double ele=1.60217653e-19;
double r=rr*1.0e-14;double nuc=nnuc*1.0e-14; // change MM to meter
double ke=(ele*ele)/(4*pai*epsi);
for (int i=1;i < 10;i++){ // repeat calculation from a to a + 10*10 (MM)
double a=aa*1.0e-14;
double rra=Math.sqrt(a*a+r*r); // rra =distance between electron and n1
double rrb=Math.sqrt(r*r+(nuc-a)*(nuc-a)); // rrb= distance between electron and n2
double poten=ke*(-1.0/rra-1.0/rrb+1.0/nuc); // poten=potential energy (J)
double ppot=poten*6.241509*1.0e18; // ppot = potential energy (eV)
double kinetic=-0.5*poten; // kinetic energy (J)
double velo=Math.sqrt((2.0*kinetic)/me); // velo = electron's velocity (m/s)
double binding=-ppot*0.5-13.606; // binding energy (eV)
double fyna=(ke*r)/(rra*rra*rra); // n1's force acting on electron in -y direction
double fxna=(ke*a)/(rra*rra*rra); // n1's force acting on electron in -x direction
double fynb=(ke*r)/(rrb*rrb*rrb); // n2's force acting on electron in -y direction
double fxnb=(ke*(nuc-a))/(rrb*rrb*rrb); // n2's force acting on electron in x direction
double internuc=ke/(nuc*nuc); // force between two nuclei
double ntwoforce = fxnb/internuc; // ratio of forces acting on n2
double debrogli = h/(me*velo); // de Broglie wavelength of electron
double wave=(2.0*pai*rra)/debrogli; // de Broglie's waves around nucleus 1
double tim=(2*pai*rra)/velo; // tim = time period of electron around n1
double vp=Math.sqrt((fynb*r)/me); // vp = electron precession veleocity around x axis
double debropre = h/(me * vp); // de Broglie wavelength of precession
double wavepre=(vp*tim)/debropre; // de Broglie's wave number of precession
double cf=(me*velo*velo)/rra; // cf= centrifugal force around n1
double fgen=ke/(rra*rra); // fgen = n1's force acting on the electron
double cosb=(rra*rra+nuc*nuc-rrb*rrb)/(2*rra*nuc); // cosi = cos b
double fnb=ke/(rrb*rrb);
fnb=(fnb*(rra-nuc*cosb))/rrb; // fnb= n2's force * cos (theta)
double fone=fnb+fgen; // total force acting on electron toward n1 (=F1)
double goux=fone*cosb; // x component of F1
double gouy=fone*Math.sqrt(1-cosb*cosb); // y component of F1
double gouxx=goux-fxnb; // subtract n2 force from goux
double fonea=Math.sqrt(gouxx*gouxx+gouy*gouy); // changed F1
double diele=Math.sqrt((nuc+a)*(nuc+a)+r*r); // between n2 and electron at position 2
double fxnbb=(ke*(nuc+a))/(diele*diele*diele); // x component of force from n2
double gouxxx=goux+fxnbb; // add n2 force to goux
double foneb=Math.sqrt(gouxxx*gouxxx+gouy*gouy); // changed F1 at position 2
fonea=(fonea+foneb)/2.0; // changed average F1
double rah=fone/cf; // rah = ratio of F1 to centrifugal force around n1
double newgoux = Math.sqrt(fonea*fonea-gouy*gouy);
double newx=(r*newgoux)/gouy;
double newradius=Math.sqrt(r*r+newx*newx); // radius influened by n2 force in x direction
double newvelo= Math.sqrt((fonea*newradius)/me);
double newtim=(2.0*pai*newradius)/newvelo; // time period
double newprewav=(vp*newtim*me*vp)/h; // new de Broglie's wave of precession
double newwav=(2.0*pai*newradius*me*newvelo)/h; // new de Broglie's waves number
System.out.printf("a:%.1f", aa);
System.out.printf(" binding energy (eV): %.4f ", binding);
System.out.printf("ratio of F1 to centrifugal: %.4f ", rah);
System.out.printf("ratio of forces acting on n2: %.3f \n", ntwoforce);
System.out.printf("waves around n1: %.3f ", wave);
System.out.printf("+ waves around x axis: %.3f ", wavepre);
int la=(int)(1000.0*wave); int lb=(int)(1000.0*wavepre);la=la+lb;
double twaves=la/1000.0;
System.out.printf(" = total waves: %.3f \n ", twaves);
System.out.printf("new waves around n1: %.3f ", newwav);
System.out.printf("+ new precession wave: %.3f ", newprewav);
la=(int)(1000.0*newwav); lb=(int)(1000.0*newprewav);la=la+lb;
twaves=la/1000.0;
System.out.printf(" = new total waves: %.3f \n ", twaves);
System.out.printf(" \n");
aa=aa+10;}}}