Why is this scipy.root code not converging?

I'm running a test problem to set up larger problems. Solving the simple unsteady heat equation via finite differences:

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from scipy.optimize import fsolve
from scipy import optimize
import pandas as pd
from numba import njit
t_max = 25
x_max = 1
N_points_1 = 61
N_points_2 = 11
dt = t_max/(N_points_1-1) 
dx = x_max/(N_points_2-1)
N_variables = N_points_1 * N_points_2
alpha = 0.01 
T_x_0 = 60 
T_x_L = 25
T_t_0 = 25
@njit
def EDP (v):
T_N =  np.ones((N_points_1, N_points_2))
#Variables
k=0
for i in range(N_points_1):
for j in range(N_points_2): 
T_N[i,j]=v[k]
k=k+1
#Model
sys = []
##Boundary Condition x = 0
for i in range(N_points_1):
BC1 = T_N[i,0] - T_x_0
sys.append(BC1)
##Boundary Condition x = L
for i in range(N_points_1):
BC2 = T_N[i,-1] - T_x_L
sys.append(BC2)
##Initial Condition
for j in range(1,N_points_2-1):
IC1 = T_N[0,j] - T_t_0
sys.append(IC1)
## Energy Balance 
for i in range(N_points_1-1):
for j in range (1,N_points_2 -1):
EB = alpha*(T_N[i,j+1]-2*T_N[i,j]+T_N[i,j-1])/dx**2 - (T_N[i+1,j]-T_N[i,j])/dt  
sys.append(EB)
sys=np.array(sys)
return sys
initial_guess = np.linspace(25,60,N_variables)        
function_compile = EDP(initial_guess)

I tested the following scipy.root algorithms:

Solution_T = fsolve(EDP,initial_guess)
Solution_T = optimize.anderson(EDP,initial_guess, verbose = True,maxiter = 1000)  
Solution_T = optimize.newton_krylov(EDP,initial_guess,verbose=True,maxiter = 1000)
Solution_T = optimize.broyden1(EDP,initial_guess,verbose=True,maxiter = 1000)
Solution_T = optimize.broyden2(EDP,initial_guess,verbose=True,maxiter = 1000)

And ONLY fsolve converges(instantly, in fact) but since i am just running a test for a larger run(which will contain about 50000 variables), fsolve won't be able to handle that and the scipy docs recommends the algorithms mentioned above, but none of them converge.

They either diverge with the objective function overflowing or remain stuck in the same value, not converging. As for why i am not using linear solvers, its because i intend to input more complex non-linear equations later.

Thanks

When it comes to solve numerically equations like this, it is usually better to express it as min(||F(x)||) instead of trying to find F(x)=0 these two formulations are equivalent when a solution exist, but if it is not possible to get an exact solution this tell you what is the closest you can get from a solution.

least_squares works for your example, and it has a few parameters that you can tune to somewhat improve the method performance

sol = optimize.least_squares(EDP, initial_guess)
print(sum(EDP(sol.x)**2))

1.383176178225409e-26