Improve compass heading calculations using 9DoF IMU

Question

I've been implementing a compass using raw data from a ICM-20948 IMU. It provides me with Accelerometer, Gyro, and Magnetometer data, but I've been using only Accelerometer and Magnetometer data to calculate heading like this:

float computeHeading(double mag_x, double mag_y, double mag_z, double accel_x, double accel_y, double accel_z)
{
    const vec3 vector_mag = create_vec3(mag_x, mag_y, mag_z);
    const vec3 vector_down = create_vec3(accel_x, accel_y, accel_z);
    const float scale = dot_vec3(vector_mag,vector_down) / dot_vec3(vector_down, vector_down);
    const vec3 vector_north = sub_vec3 (vector_mag , scale_vec3(scale, vector_down));
    return atan2(vector_north.x, vector_north.y) * 180 / M_PI;
}

It works well when doing the YAW rotation, but everything falls apart once I PITCH or ROLL the IMU for just a couple of degrees. Here is where I got the algorithm from.

Is there any way I can make the calculation more robust? (maybe using gyro)