Mira si las fórmulas son correctas. Obtengo valores extraños y, al menos la fórmula de la ecuación de 2º grado que has dado en el enunciado está mal.

Estoy oxidado en matemáticas.

Tengo los siguientes resultados:
d(20º) = 1837.556662
d(40º) = 4492.157928
d(60º) = 5087.585718
d(80º) = 2250.260680

¿Tienes los resultados de la ecuación para compararlos?

Replanteé la fórmula en base a la física básica de instituto. Hace tanto que no toco todo esto que no puedo asegurar que esté bien, por eso te pedía los resultados, si los tenías.

Partiendo de la formula de la aceleración: a = (vf - vo) / (tf - to)
y de la fórmula del espacio recorrido por una velocidad constante: e = v * t

donde
a : aceleración en m/s^2
v : velocidad en m/s
t : tiempo en s

Siendo:
vx : velocidad horizontal
vy : velocidad vertical
vo : velocidad inical
A : ángulo de disparo
ym : altura del cañón respeto al suelo
g : aceleración de la gravedad = 9.8m/s^2

vx = vo * cos(A)
vy = vo * sin(A)

1. ALTURA MÁXIMA
g = (vf - vo) / (tf - to) = -vo / tf

La aceleración es negativa puesto que va frenando. Consideramos pues -g
-g = -vo / tf
tf = -vo / -g = vo / g   =>   t1 = vy / g

e = v * t    =>    h = vy * t1      ( h es la altura máxima que adquiere el proyectil respecto al cañón )
h = vy * (vy / g) = vy^2 / g

ht = ym + h    (ht es la altura máxima con respecto al suelo)

2. TIEMPO DE CAÍDA
g = (vf - vo) / (tf - to)
v = e / t

g = (ht / tp) - 0 / 1 - 0 = ht / tp    ( tp es el tiempo desde la máxima altura hasta que toca suelo )
tp = ht / g

3. ESPACIO HORIZONTAL RECORRIDO
t = t1 + tp    ( t es el tiempo de vuelo del proyectil: el tiempo de subida más el tiempo de bajada )
e = vx * t

4. SUSTITUYENDO
e = vx * t
e = vx * (t1 + tp)
e = vx * (vy/g + ht/g) = vx * [(vy + ht) / g] = vx * (vy + ht) / g
e = vx * (vy + ym + h) / g
e = vx * (vy + ym + vy^2 / g) / g

De esta forma el programa resulta en lo siguiente:

#include <stdio.h>
#include <math.h>
 
#define M_PI 3.14159265359
#define DEG2RAD(X) (X)*M_PI/180
 
int main() {
    const double g = 9.8;
    const double Vo = 105;
    const double Ym = 15;
    double Vx;
    double Vy;
    double angulo;
    double angulof;
    double d;
 
    for(angulo = 20; angulo <= 80; angulo += 20) {
        angulof = DEG2RAD(angulo);
        Vx = Vo * cos(angulof);
        Vy = Vo * sin(angulof);
        d = Vx * (Vy + Ym + Vy*Vy/g) / g;
        printf("d(%02.0fº) = %f\n", angulo, d);
    }
    return 0;
}