feedback: solve_bvp boundary conditions

Question

This is my first time posting here so pardon my lack of experience formatting these posts.

Basically I have a system of ODEs I'm solving with solve_bvp and my results are good, but not great. I was wondering if it's an issue with my boundary conditions. Unfortunately, for the functions of Tg I only have values for Tg2(a) and Tg2(b), and nothing for Tg1(a), except that Tg2 is the derivative of Tg1 (resulting from making the second order Tg into 2 first order ODEs: Tg1 and Tg2).

Should I rest easy knowing this is as good as it gets with the BCs I have?

#System of ODEs
def bvp(x, T): #apparently the bvp solver only works with x, so let z=x for the purpose of this code
    Ts1 = T[0]
    Ts2 = T[1]
    Tg1 = T[2]
    Tg2 = T[3]
    a = 1 #simplified for the purpose of posting here

    dTs1dx = Ts2
    dTs2dx = Ts2*a + Ts1*a + a + a*(Ts1-Tg1) + a*(Ts1-a)
    dTg1dx = Tg2
    dTg2dx = Tg2*a + Tg1*a - a*(Ts1-Tg1) + a(Tg1-a)
    return [dTs1dx, dTs2dx, dTg1dx, dTg2dx]

#Boundary Conditions
def bcs(Ta, Tb): #Boundary Residuals of bvp at boundaries a and b, where a=0 and b=Lr
    Ts1a, Ts2a, Tg1a, Tg2a = Ta #Ts, dTs, Tg, and dTg values at boundary a=0
    Ts1b, Ts2b, Tg1b, Tg2b = Tb  #T and dT values at boundary b=Lr
    return [Ts1a - (25+273), Ts2b - (6530.34), Tg2a, Tg2b - (5450.2)]