Double scalars

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
import copy

 
np.seterr(all='raise')

 
subdivision = 1000
a = 0.3
sigma = 5.67*(10**(-8))
phi0 = 1360
epsilon = 0.25
l = 40
dtheta = (np.pi)/(2*subdivision)
 

angles = np.linspace(0, np.pi/2, num=subdivision)
 

def phi(y0, theta) :
    x, xp = y0
    xpp = (np.tan(theta))*xp + \
        ((x**4)*sigma*(1+epsilon)/(l*(np.cos(theta)))) \
            - ((1-a)*(1-epsilon)*phi0/((np.pi)*l))
    return(xp, xpp)

 
x0 = 229.11

x = odeint(phi, [x0,0], angles)

T = (x[::,0]**4 + ((np.cos(angles))*phi0*(1-a)/((np.pi)*sigma)))**(1/4)

 

def A (temp) :
    if temp<273.15 :
        return (0.8)
    else :
        return (0.3)

 

def phi2(y0, theta) :
    i = int(theta/(np.pi/(2*subdivision)))
    x, xp = y0
    xpp = (np.tan(theta))*xp + \
        ((x**4)*sigma*(1+epsilon)/(l*(np.cos(theta)))) \
            - ((1-A(T[i]))*(1-epsilon)*phi0/((np.pi)*l))
    return(xp, xpp)
 

Ta0 = odeint(phi2, [x0,0], angles)
 
T0 = (Ta0[::,0]**4 + ((np.cos(angles))*phi0*(1-a)/((np.pi)*sigma)))**(1/4)

def lissage (T) :
    j = 1
    for i in range (1, len(T)-1) :
        d2T = (T[i+1]-2*T[i]+T[i-1])/((dtheta)**2)
        if j==1 and d2T >= 0 :     
            j = i
    deltaT = T[j] - T[j+1]
    for i in range (j, len(T)) :
        T[i] = T[i-1] - deltaT
    return(T,j)


(T,j) = lissage(T0)
Ta = np.copy(T)
for i in range(0, subdivision):
    Ta[i]=((T[i]**4)-np.cos(angles[i])*phi0*(1-a)/((np.pi)*sigma))**(1/4)
    
    
def evolution (T0, Ta, iterations) :
    for k in range (1, iterations +1) :
        T1 = np.copy(T0)
        Ta1 = np.copy(Ta)
        for i in range (1, subdivision-1) :
            dTa = (Ta[i+1]-Ta[i])/dtheta
            d2Ta = (Ta[i+1]-2*Ta[i]+Ta[i-1])/((dtheta)**2)
            Ta1[i] = (((l*((np.cos(angles[i]))*d2Ta - (np.sin(angles[i]))*dTa)+2*(1-A(T0[i]))*phi0*np.cos(angles[i])/np.pi)/(sigma*(1+epsilon))) - (1-A(T0[i]))*phi0*np.cos(angles[i])/(sigma*np.pi))**(1/4)
            T1[i] = ((Ta1[i]**4)+((1-A(T0[i]))*phi0*np.cos(angles[i])/(sigma*np.pi)))**(1/4)
        T0 = T1
        Ta = Ta1
        plt.plot(angles, T0)
        plt.show()

evolution(T, Ta, 100)

((l*((np.cos(angles[i]))*d2Ta - (np.sin(angles[i]))*dTa)+2*(1-A(T0[i]))*phi0*np.cos(angles[i])/np.pi)/(sigma*(1+epsilon))) - (1-A(T0[i]))*phi0*np.cos(angles[i])/(sigma*np.pi)