% plot of intensity 
T1=6000; T2=4000;
h=6.625*10^(-34);
k=1.38*10^(-23);
c=3*10^8;
lambda = .03*(1:100)*10^(-6);
int1 = pi*(2*h*c^2)*lambda.^(-5)./(exp(h*c./(k*T1*lambda)) - 1);
int2 = pi*(2*h*c^2)*lambda.^(-5)./(exp(h*c./(k*T2*lambda)) - 1);
plot(lambda,int1,'.',lambda,int2,'.');