solve heat ode

from scipy.integrate import solve_bvp
from numpy import *
import matplotlib.pyplot as plt 
####################################
# 1/r (k t' + kr t") + e' = 0
# Sea y0 = t,  y1 = t'
# 
# [ y0 ]'  =  [ y1               ]
# [ y1 ]      [( -e' r - ky0)/kr ]
####################################

## y = temp, r = radio
## y = y(r)

r_saltos = 100 ## subdivisiones de radio
r = linspace(0,0.0025,r_saltos) ## rango del radio
y = array([r, full_like(r,2)])


def edo(r,y):   ## sistema de primer orden de edos, variable (indep, dep)
                
    e = (10 * 220) / (pi *(0.005 / 4)**2) # e'(t) = P / vol            
    k = -0.076 * y[0] + 450.158
    
    return vstack([ y[1] , (-1*e*r - k*y[0]) / (k*r) ])

def bc(ya,yb):  ### condiciones de contorno
    return array([ya[1],yb[0]-573])