Boundary Value Problem with ModellingToolkit or DiffEqOperators

At first I thought that there was a bug in the packages but probably I am doing something (completely wrong). I am trying to run the bellow spherical coordinates system with Boundary Values. But not sure if it even makes sense…

@parameters r
@variables T(..)
D = Differential(r)

kc = 25.
ks = 385.
Tf = 350.
Rc = 0.005
Rs = 1.5 * Rc
h = 100.
Sc = 0.
Tinit = 298.
	
# 1D ODE and boundary conditions
eq = [0 ~ (1/r^2) * D( (1/r^2) * kc * D( T(r) ) ) + Sc,
      0~ D(T(0)),
      0~ D(T(Rs)) - h * (T(Rs) - Tf)]  

# ODE system
sys = ODESystem(eq, T)
p=[]
T0 = [T => Tinit]
rspan=(0.0 , Rs)
prob = ODEProblem(sys, T0, rspan, p)

# Solve ODE problem
sol = solve(prob)

BVPs are not supported in MTK right now. Please open an issue. You can use BVProblem from DifferentialEquations.jl directly for now though.

https://diffeq.sciml.ai/stable/tutorials/bvp_example/

thanks for the update. For reference here is the issue
https://github.com/SciML/ModelingToolkit.jl/issues/924

BTW, is there a way to set BCs for the derivatives (Neumann, Robin)?

Yes, in PDESystem form.