restart; with(plots):

Earth's heat balance

define system...

f := T -> -epsilon*(T+273)^4+w(T)*cos(theta);

w := T -> a+b*erf(c*T);

df := simplify(D(f));

a := 4.4e9; b := 1.3e9; c := .2; epsilon := .5;

plot([w(T)], T = -30 .. 30, labels = (["T", "w(T)"]));

determine fixed points

theta := 0;

plot([f(T)], T = -60 .. 60, labels = (["T", "dT/dt"]));

only one fixed point

fp[1] := fsolve(f(T) = 0, T = 50);

sigma[1] := evalf(df(fp[1]));

So T = 53.7C is stable

theta := 1/4*Pi;

plot([f(T)], T = -60 .. 60, labels = (["T", "dT/dt"]));

three fixed points

fp[1] := fsolve(f(T) = 0, T = -50);

fp[2] := fsolve(f(T) = 0, T = 0);

fp[3] := fsolve(f(T) = 0, T = 50);

sigma[1] := evalf(df(fp[1]));

sigma[2] := evalf(df(fp[2]));

sigma[3] := evalf(df(fp[3]));

So: T = -15.6C is stable, T=-2.1C unstable, T = 26.6C is stable

The domains of attraction are:

fp[1]: -273C <= T < fp[2]

fp[2]: empty

fp[3]: fp[2] < T < infinity

theta := 1/3*Pi;

plot([f(T)], T = -60 .. 60, labels = (["T", "dT/dt"]));

only one fixed point

fp[1] := fsolve(f(T) = 0, T = 50);

sigma[1] := evalf(df(fp[1]));

So T=-37 C is stable

Potential Function: we determine V(x) such that dV/dx=-f(x):

theta := 1/4*Pi;

V := unapply(-int(f(T), T), T);

plot([V(T)], T = -60 .. 60, labels = (["T", "V(T)"]));

Evolution in time for two initial conditions

epsilon := 'epsilon'; theta := 'theta';

eq := diff(T(t), t) = f(T(t));

epsilon := .5; theta := 1/4*Pi;

ic1 := T(0) = 100;

ic2 := T(0) = -100;

sol1 := dsolve({ic1, eq}, T(t), type = numeric);

sol2 := dsolve({eq, ic2}, T(t), type = numeric);

tinf := 1e-6;

pl1 := odeplot(sol1, [t, T(t)], 0 .. tinf, color = red, legend = "hot");

pl2 := odeplot(sol2, [t, T(t)], 0 .. tinf, color = blue, legend = "cold");

display(pl1, pl2, labels = (["t", "T[C]"]));

equilibrium solutions as a function theta

thetamin := 0.0;

thetamax := Pi/2;

n := 50;

tinf := 1e-6;

loop...

bif1 := NULL: bif2 := NULL:

for i from 0 to n do

theta := thetamin + i*(thetamax-thetamin)/n:

sol1 := dsolve({eq,ic1},T(t), type=numeric);

sol2 := dsolve({eq,ic2},T(t), type=numeric);

T1 := subs(sol1(tinf),T(t));

T2 := subs(sol2(tinf),T(t));

bif1 := bif1, [theta/Pi*180, T1]:

bif2 := bif2, [theta/Pi*180, T2]:

od:

plot...

pl1 := plot([bif1],labels=["theta","T[C]"],legend="hot branch",color=red):

pl2 := plot([bif2],labels=["theta","T[C]"],legend="cold branch",color=blue):

display(pl1,pl2);

To get an indication of the critical theta values we look at the graphs of f(T)

theta := 'theta';

plot([subs(theta = (38*Pi)*(1/180), f(T)), subs(theta = (57*Pi)*(1/180), f(T))], T = -40 .. 40,legend = (["theta=38", "theta=57"]));

soltheta1 := fsolve({f(T)=0,df(T)=0},{T=-5,theta=38*Pi/180});

soltheta2 := fsolve({f(T)=0,df(T)=0},{T=8,theta=57*Pi/180});

thetacrit[1] := evalf(180*subs(soltheta1, theta)/Pi);

thetacrit[2] := evalf(180*subs(soltheta2, theta)/Pi);