double dt = 0.1;
double  t = 0.0;

Integrator integrator;

double[] y;
double[] y0;
double[] y_alt;

void setup()
{
    size(800,600);  
    
    Modell modell = new Modell();
    y = new double[modell.anzahl];
    y0 = new double[] {1.0,0.0};
    y[0] = y0[0];
    y[1] = y0[1];
    y_alt = new double[modell.anzahl];
    integrator = new Integrator(modell);

    frameRate((int)(1.0/dt));
}

void draw()
{
    clear();
    fill(255,0,0);
    ellipse(width/2,height/2+(float)(100.0*y[0]),20,20);


    for(int i=0;i<y.length;i++)
        y_alt[i] = y[i];
  
    integrator.zeitschritt(y,y_alt,t,dt); //Ergebnis wird in y gelegt
  
  
    t+=dt;   
}

Code 0-1: simulation.pde

public class Modell
{
    double D = 0.1;
    double C = 1.0;
    double m = 1.0;
  
    public int anzahl = 2; //Anzahl der Gleichungen , ANPASSEN!
  
    void steigung(double[] f, double[] y,double t)
    {
        double x = y[0];
        double v = y[1];
 
        f[0] = v;
        f[1] = -(C/m)*x -(D/m)*v; 
    }  
}

Code 0-2: Modell.pde

public class Integrator
{
    int anzahl;
  
    double[] k1,k2,k3,k4;
    
    double[] yhilf;
  
    Modell modell;
  
    public Integrator(Modell modell)
    {
        this.modell = modell;
        anzahl = modell.anzahl;
        k1 = new double[anzahl];      
        k2 = new double[anzahl];      
        k3 = new double[anzahl];      
        k4 = new double[anzahl];
  
        yhilf = new double[anzahl];            
    }  
    
    public void zeitschritt(double[] y,double[] y_alt, double t, double dt)
    {
        modell.steigung(k1,y_alt,t);
        for(int i=0;i<anzahl;i++)
            yhilf[i] = y_alt[i] + 0.5*dt*k1[i];
            
        modell.steigung(k2,yhilf,t+0.5*dt);            
        for(int i=0;i<anzahl;i++)
            yhilf[i] = y_alt[i] + 0.5*dt*k2[i];

        modell.steigung(k3,yhilf,t+0.5*dt);            
        for(int i=0;i<anzahl;i++)
            yhilf[i] = y_alt[i] + dt*k3[i];

        modell.steigung(k4,yhilf,t+dt);

        for(int i=0;i<anzahl;i++)
            y[i] = y_alt[i] + (dt/6.0)*(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i]);        
     }
}

Code 0-3: Integrator.pde