I have kinda solved the problem - instead of using quaternions by themselves I have calculated a new vector between the forward vector and the vector pointing to the model using a modified spherical interpolation method (instead of passing a fraction I pass an angular speed) and then calculated a lookAt matrix like in this article http://stackoverflow.com/questions/349050/calculating-a-lookat-matrix . Now I'm just still wondering what is wrong with the quaternion method. I hope that my actual solution doesn't introduce Gimbal Lock problem.

Code of my current rotateTowards function:

	rotateTowards: function(target, angularSpeed) {
		if (this.position.isEqual(target)) { // If is on the same position
			return;
		}

		var vecStart = this.rotation.getForwardVector();
		var vecTarget = this.position.sub(target).normalize();
		var angleLeft = Math.acos(vecStart.dot(vecTarget));
		var fullAngle = angularSpeed * deltaT;

		if (Math.areScalarsSimilar(angleLeft, 0)) { // If already rotated towards the target
			return;
		}
		else if (Math.areScalarsSimilar(angleLeft, Math.PI)) { // If vectors are opposite
			// Turn around up-axis a bit
			this.rotation.multiply(new Quaternion().axisToQuaternion(fullAngle, upVec()));
			MVMatrix.setRotation(this.rotation.quaternionToMatrix());
		}
		else { // If vectors casually differs
			if (angleLeft > fullAngle) { // If able to make another full rotation step
				var zAxis = new vec3().angularSlerp(vecStart, vecTarget, fullAngle);
				var xAxis = upVec().cross(zAxis).normalize();
				var yAxis = zAxis.cross(xAxis);
				MVMatrix.make(
					xAxis.x,			xAxis.y,			xAxis.z,
					yAxis.x,			yAxis.y,			yAxis.z,
					zAxis.x,			zAxis.y,			zAxis.z,
					this.position.x,	this.position.y,	this.position.z);

				this.rotation.matrixToQuaternion(MVMatrix);
				this.pitch = this.rotation.getPitch();
				this.yaw = this.rotation.getYaw();
				this.adjustAngles();
			}
			else { // If the angle left is smaller than the angular speed
				// Set the final rotation
				this.lookAt(target);
			}
		}
	},

And the simple angularSlerp function:

	angularSlerp: function(v0, v1, angularSpeed) {
		var angleBetween = Math.acos(v0.dot(v1));

		v0.multiplyByScalar(Math.sin(angleBetween - angularSpeed)/Math.sin(angleBetween));
		v1.multiplyByScalar(Math.sin(angularSpeed)/Math.sin(angleBetween));
		return v0.add(v1);
	},