function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity)
    if this.Gravity
        g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
    else
        g = 0;
    end

    % Probability of Background Collisions
    collision = rand(1);
    CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime);

    if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision

        ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9);

        for i=1:int64(this.SimulationTime/this.TimeStep)
            
            ParticleTrajectory(i,1:3) = InitialPosition;
            ParticleTrajectory(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;
            ParticleTrajectory(i,7:9) = (ga1+2*(ga2+ga3)+ ga4)./6;
        end

        FinalDynamicalQuantities = ParticleTrajectory(end,:);

    else
        ParticleTrajectory       = zeros(int64(this.SimulationTime/this.TimeStep),9);
        FinalDynamicalQuantities = -500*ones(1,9); % -500 is a flag for giving up this trajectory
    end
end