Scipy.root not converging even when provided with initial guesses very close to solution

Question

I've made a previous question here and also in SO wondering why only the fsolve solver converges for the simple one dimensional unsteady conduction problem

$$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2T}{\partial x^2}$$

$\alpha = 0.001$

$T(x=0,t) = 60°C$

$T(x=L,t) = 25°C$

$T(x,t=0) = 25°C$

One of the answers alluded to initial values, from which i started to test and, to make things easier, i changed the value of alpha to

$\alpha = 0.000001$

Which means that the result profile shall be very close to simply

$T(x=0,t) = 60°C$

$T(x=L,t) = 25°C$

$T(x,t) = 25°C$

Which is what i used as initial condition for this new code

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

t_max = 25
x_max = 1
N_points_1 = 100
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.000001 
T_x_0 = 60 
T_x_L = 25
T_t_0 = 25


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  #partial(T_N,t_1,1)[i,j] - alpha*partial_2(T_N,x,2)[i,j] 
            sys.append(EB)
    
    sys=np.array(sys)

    
    return sys

initial_guess = 25*np.ones((N_points_1,N_points_2))
initial_guess[:,0] = T_x_0
initial_guess[:,-1] = T_x_L
initial_guess = initial_guess.flatten()
print(np.sum(EDP(initial_guess)))
0.3464999999999999

As you can see the value of the sum of the values of each equation is already close to 0

fsolve does converge, but even so, it takes a long number of iterations. Others either do not converge or straigth up diverge

Solution_T = fsolve(EDP,initial_guess,maxfev = 1200)

Solution_T2 = optimize.anderson(EDP,initial_guess, verbose = True,maxiter = 1000)
104:  |F(x)| = 2.36365e+93; step 1
OverflowError: (34, 'Result too large')

Solution_T3 = optimize.newton_krylov(EDP,initial_guess,verbose=True,maxiter = 1000)
999:  |F(x)| = 0.00372199; step 1
NoConvergence: [60.         24.99999991 24.99999972 ... 24.99990242 24.99977307
 24.99968681]

Solution_T4 = optimize.broyden1(EDP,initial_guess,verbose=True,maxiter = 1000)
659:  |F(x)| = 4.37351e+151; step 1
OverflowError: (34, 'Result too large')

Solution_T5 = optimize.broyden2(EDP,initial_guess,verbose=True,maxiter = 1000)
999:  |F(x)| = 2.81888e+10; step 1.74723e-09
NoConvergence: [ 6.00000000e+01 -2.63924280e+06 -2.98161396e+04 ...  3.39690111e+08
  2.14121070e+07 -2.54528649e+08]

So, what actually is happening here? I can only converge these algorithms when i feed the solution of fsolve into them, from which the objective function is already below the tolerance for convergence. I know that fsolve did converge, but i am just running tests for much larger system of equations, from which the large scale solvers, those above besides fsolve, are required.