One problem is a typo in ‘E2’, where there’s an ‘e’ wanting in ‘tprima’. With that and creating ‘E1E2’ from both functions:
alphaS=0.11; ta=0.53; alphaA=0.3; tprimea=0.06; alphaprimea=0.31;
C=2.5;
So=6.28e7; %Solar Constant, p.110
sigma=5.67e-8; %Boltzman's Constant, p. 110
E1=@(Ts,Ta)ta*(1-alphaS)*So/4-C*(Ts-Ta)-sigma*Ts^4*(1-alphaprimea)+...
sigma*Ta^4;
E2=@(Ts,Ta)(1-alphaA-ta*(1-alphaS))*So/4+C*(Ts-Ta)+sigma*Ts^4*...
(1-tprimea-alphaprimea)-2*sigma*Ta^4;
E1E2 = @(B) [E1(B(1),B(2)); E2(B(1),B(2))];
B0 = [290; 250];
TsTa = fsolve(E1E2, B0)
TsTa =
4.5604e+003
3.5995e+003
Since this may not be the result you want, you may need to experiment with it, particularly with respect to a different ‘B0’. (This code assumes both equations equate to zero, as fsolve requires.)
