class Simulation{

private DrawSim plot;
private Quantum oneD;
public double energy;
public boolean parity;

Simulation(DrawSim plotting, Quantum oneDlocal){
oneD=oneDlocal;
plot=plotting;
parity=true;
}

public void go(){
int i;
double y[] = new double[2];
double yout[]= new double[2];
double x,dx;
int nsteps;

dx=0.01;
nsteps=(int)(5.0/dx);

plot.clearScreen();

// initial conditions

if(parity){
y[0]=1.0;
y[1]=0.0;}
else{
y[0]=0.0;
y[1]=2.0;}

x=0.0;

// integrate 

                for(i=0;i<=nsteps;i++){
                        rk4(y,yout,x,dx);
                        plot.addToScreen(x+dx,yout[0]);
                        y[0]=yout[0];
                        y[1]=yout[1];
                        x=x+dx;
                }

return;
}


private void f(double x, double y[], double dy[]){
        int i;

        for(i=0;i<2;i++){
                if(i==0) dy[0]=f0(x,y);
                if(i==1) dy[1]=f1(x,y);
        }

        return;
}

// derivative of psi 

private double f0(double x, double y[]){
        double dum;

        dum=y[1];
        return(dum);
}

private double f1(double x, double y[]){
        double dum;

        dum=(0.25*x*x-energy)*y[0];
        return(dum);
}

private void rk4(double y[], double yout[], double x, double dx){
        double k1[],k2[],k3[],k4[];
        double yt[],dy[];
        int i;
 
        k1 = new double[2];
        k2 = new double[2];
        k3 = new double[2];
        k4 = new double[2];
        yt = new double[2];
        dy = new double[2];

        f(x,y,dy);
        for(i=0;i<2;i++){
                k1[i]=dx*dy[i];
                yt[i]=y[i]+0.5*k1[i];
        }

        f(x+dx*0.5,yt,dy);
        for(i=0;i<2;i++){
                k2[i]=dx*dy[i];
                yt[i]=y[i]+0.5*k2[i];
        }

        f(x+dx*0.5,yt,dy);
        for(i=0;i<2;i++){
                k3[i]=dx*dy[i];
                yt[i]=y[i]+k3[i];
        }

        f(x+dx,yt,dy);
        for(i=0;i<2;i++)
                k4[i]=dx*dy[i];

        for(i=0;i<2;i++)
                yout[i]=y[i]+(1.0/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

        return;
}

}