Trying to print and plot two parameters of multivariate functions

edit

I am trying to numerically find the values of two parameters that will make three multivariate functions true. I am trying to use for loops and if statements to do this. I need this for a university project I am working on. I am completely new to Python, and the syntax of the Sage notebook I have created does not seem to be working. Below is what I typed. Also, here is the link to the notebook I created.

var('V_m,phi_0,f_0,Lambda,sigma_q,T_i,t,T_2,I_0,sigma_B,T_amb,eta')
V_m=3.06*10^-4
T_i=1300
Lambda=4.5
phi_0=0.5
f_0=srange(0.1,100,step=0.1,universe=None)
T_2=srange(500,2000,step=0.1,universe=None)
sigma_q=srange(0.5*10^12,3.5*10^12,step=10^6,universe=None)
I_0=700000
T_amb=300
pi_n=pi.numerical_approx()
sigma_B=5.67*10^-8
for i in f_0:
    for j in T_2:
        for k in sigma_q:
            R_0=n(((V_m/((1-phi_0)*float(pi)*i))^(1/3))*sqrt(phi_0))
            H_0=n(((V_m/((1-phi_0)*float(pi)*i))^(1/3)))
            L_0=n((((V_m*(i)^2)/((1-phi_0)*float(pi))))^(1/3))
            T_R0L0=n(T_i+((k)/(2*Lambda))*(1/2)*(((R_0)^2)-((H_0)^2))+((H_0)^2)*ln((H_0)/(R_0)))
            T_MAXR0=n((((I_0)/(float(pi)*((H_0)^2)*sigma_B))+((T_amb)^4)-((k)/(sigma_B))*((((H_0)^2)-((R_0)^2))
            /R_0))^(1/4))
            T_Mean=n((T_MAXR0+T_R0L0)/2)
            T_R00=n(((1/2)*(((I_0)/(float(pi)*((H_0)^2)*sigma_B))+((T_amb)^4)+((T_Mean)^4)-((k)/(sigma_B))*
            ((((H_0)^2)-((R_0)^2))/R_0)))^(1/4))
            A1=n((-1/2*I_0*(i - sqrt(i^2 + 4*phi_0))*i))
            A2=n(V_m * k + 1/30*(15*j^4*i^6 + 15*T_amb^4*i^6 - 6*(i^6 - 20*(3*i^2*phi_0^2
            - 4*phi_0^3)*arcsinh(1/2*i/sqrt(phi_0)))*T_R00^4 - 6*(i^6 + 40*(3*i^2*phi_0^2 -
            8*phi_0^3)*arcsinh(1/2*i/sqrt(phi_0)))*T_R00^3*T_R0L0 - 6*(i^6 -120*(i^2*phi_0^2 -
            4*phi_0^3)*arcsinh(1/2*i/sqrt(phi_0)))*T_R00^2*T_R0L0^2 - 6*(i^6 + 40*(3*i^2*phi_0^2 -
            8*phi_0^3)*arcsinh(1/2*i/sqrt(phi_0)))*T_R00*T_R0L0^3 - 6*(i^6 - 20*(3*i^2*phi_0^2 - 
            4*phi_0^3)*arcsinh(1/2*i/sqrt(phi_0)))*T_R0L0^4 - 3*(5*j^4*i^5 + 5*T_amb^4*i^5 - 2*(i^5 + 
            3*i^3*phi_0 - 44*i*phi_0^2)*T_R00^4 - 2*(i^5 - 2*i^3*phi_0 + 176*i*phi_0^2)*T_R00^3*T_R0L0 - 
            2*(i^5 - 2*i^3*phi_0 - 264*i*phi_0^2)*T_R00^2*T_R0L0^2 - 2*(i^5 - 2*i^3*phi_0 +
            176*i*phi_0^2)*T_R00*T_R0L0^3 - 2*(i^5 + 3*i^3*phi_0 - 44*i*phi_0^2)*T_R0L0^4)*\
            sqrt(i^2 + 4*phi_0) + 32*(144*T_R00^2*T_R0L0^2*i*phi_0^2 - (5*i^3*phi_0 - 24*i*phi_0^2)\
            *T_R00^4 + (5*i^3*phi_0 - 96*i*phi_0^2 ...