sry i just wake up.

We suppose an angle w, the angle between the X direction and the direction of the vector V0

So we got Vx0=V0*cos(w) and Vy0=V0*sin(w)

We calculate the differance of position between the point P(X,Y) at instant t on the parabola and the initial point P0(X0,Y0) at instant t0

We obtain: X-X0=Vx0*t+X0-Vx0*t0-X0 but t0=0 so (X-X0)/(Vx0)=t [equation 1]

We calculate the same for Y-Y0

Y-Y0=(1/2)*g*t²+Vy0*t+Y0-(1/2)*g*t0²-Vy0*t0-Y0 but t0=0

so Y=(1/2)*g*t²+Vy0*t+Y0

You replace t by his value in the [equation 1] and after some simplification you got:

Y=(1/2)*g*(X-X0)²/(Vx0)²+Vy0*(X-X0)/(Vx0)+Y0

If you consider the origin at P0(X0=0,Y0=0) it will simplify the calculation after(if not use the same technique as below but with a bigger equation):

Y=(1/2)*g*(X)²/(Vx0)²+Vy0*(X)/(Vx0)

And with the power of trigo:

1/(Vx0)²=1/(V0*cos(w))²=(1+[tan(w)]²)/V0²

and Vy0/Vx0=V0*sin(w)/V0*cos(w)=tan(w)

we replace them in the equation:

Y=(1/2)*g*(X)²(1+[tan(w)]²)/V0²+X*tan(w)

Then you use the by point P1 the parabola pass by:

Y1=(1/2)*g*(X1)²(1+[tan(w)]²)/V0²+X1*tan(w)

And i let you do the rest, it is now a polynome and 2 degre (solve it conserdiring tan(w)=T so you have to solveY1=(1/2)*g*(X1)²(1+T²)/V0²+X1*T

Then when it solve use the value S you find and do w=arctan(S)

Caution: you will have these case:

-if P1 is out of range: polynome insolvable

-if P1 is at max range : 1 solution (noramlly it will be w=45°)

-if P1 is between max range and initial position so 2 solutions (one with w>45° and other with w<45°)

Hope it is what you want, if you prefer Vx0 and Vy0 value just do the calculation with the value of w in the first equation:Vx0=V0*cos(w) and Vy0=V0*sin(w)

Tell me if you have questions.