Y = 2021:2100;
C = zeros(1, 80);
M = zeros(1, 80);
N = zeros(1, 80);
T = zeros(1, 80);


for i = 1:length(Y)
    C(i) = 424 + (Y(i) - 2021)*2.4;
    M(i) = 1911.8 + (Y(i) - 2021)*12;
    N(i) = 337.5 + (Y(i) - 2021)*1;

    F1 = 5.35*log(C(i)/278);
    F2 = 0.036*(sqrt(M(i)) - sqrt(700));
    F3 = 0.12*(sqrt(N(i)) - sqrt(270));

    T(i) = ((4 * (F1+F2+F3) + 1361.2 * (1-0.31))/(4 * 0.61 * 5.67e-8))^0.25 - 273.15;

endfor

plot(Y, T,'ko');
xlabel('Year', 'FontSize', 16)
ylabel('Temperature (^oC)', 'FontSize', 16)
axis('tight')
set(gca,'FontSize',16)