What this problem shows is that Green's function is a generalization of "variation of parameters".
If y1 and y2 are complimentary solutions, then you seek a specific solution of the form y(x)= u(x)y1(x)+ v(x)y2(x). Of course, there are many function, u and v, that will satisfy that. Differentiating, y'= u'y1+ uy1'+ v'y2+ vy2', by the product rule. Because there are many possible solutions, we narrow the search and simplify the problem by requiring that u'y1+ v'y2= 0. Now, we have y'= uy1'+ vy2' so y"= u'y1'+ uy1"+ v'y2'+ vy2". Putting that into the original differential equation, y"+ py'+ qy= u'y1'+ uy1"+ v'y2'+ vy2"+ puy'1+ pvy2'+ quy1+ qvy2= (uy1"+ puy1'+ quy1)+ (vy2"+ pvy2'+ qvy2)+ (u'y1'+ v'uy2')= g(x).
But uy1"+ puy1'+ quy1= u(y1"+ py1'+ qy1)= 0 because y1 satisfies the homogenous equation and vy2"+ pvy2"+ qvy2= v(y1"+ py1'+ qy1)= 0 because y2 satisfies the homogenous equation. That leave u'y1'+ v'y2'= g(x).
Together with u'y1+ v'y2= 0, that give two equations we can solve, algebraically, for u' and v'. Do that, then integrate to get the formula you want.