I am trying to implement a solution of a dynamic EDP system with two independent variables (t,z). I already discretized the problem and fed into scipy.optimize.newton_krylov. The overall system is non linear because of some extra non linear correlations with exponential functions needed to calculate parameters that feed the EDPs and also change dynamically.
The overall system contains 342108x342108 equations x variables and so far its taking 4+ hours to converge to no avail. I know the convergence is slow because i have the exact same model on an external(paid) program that solves EDPs and at the slowest it takes 40 minutes to converge, using the exact same discretization method.
In fact, the external program can easily deal with 3million equations on finer grids and solve it in 5 minutes for the fastest cases.
The model itself is enormous, it has about 100 equations(without discretization) and each differential equation is pretty big. Hence why i did not paste the code here.
in terms of initial guesses, I purposely put values that are close to the final solutions based on an external run, still takes far too long to converge. So what can i do?