function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity)
    
    if this.Gravity
        g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
    else
        g = 0;
    end
    
    ParticleDynamicalQuantities   = zeros(int64(this.SimulationTime/this.TimeStep),6);
    for i=1:int64(this.SimulationTime/this.TimeStep)

        ParticleDynamicalQuantities(i,1:3) = InitialPosition;
        ParticleDynamicalQuantities(i,4:6) = InitialVelocity;

        rt = InitialPosition;
        vt = InitialVelocity;

        ga1 = this.calculateTotalAcceleration(rt,vt) + g;
        gv1 = vt .* this.TimeStep;
        rt  = rt + 0.5 * gv1;
        vt  = vt + 0.5 * ga1 .* this.TimeStep;

        ga2 = this.calculateTotalAcceleration(rt,vt) + g;
        gv2 = vt .* this.TimeStep;
        rt  = rt + 0.5 * gv2;
        vt  = vt + 0.5 *ga2 .* this.TimeStep;

        ga3 = this.calculateTotalAcceleration(rt,vt) + g;
        gv3 = vt .* this.TimeStep;
        rt  = rt + 0.5 * gv3;
        vt  = vt + ga3 .* this.TimeStep;

        ga4 = this.calculateTotalAcceleration(rt,vt) + g;
        gv4 = vt .* this.TimeStep;

        InitialPosition  = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
        InitialVelocity  = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
        %Acceleration    = (ga1+2*(ga2+ga3)+ ga4)./6;
    end
end