The ‘Ps’ integration produceds onlly one value, the value of the integral between the limits of integration.
To get a vector for it, either use a for loop or arrayfun (which is a for loop in disguise) as I did here.
% clc
% clear all
Fs= 2.16*10^-6*pi; % Geometrical Factor
h= 4.136*10^-15; % Plancks Constant
c= 3*10^8; % Speed of light
K = 8.61*10^-5; % Boltzmanns Constant
Ts=5760; % Temparature of the sun
% E=0:0.0001:4;
E = linspace(0, 6, 250);
bs=((2*Fs)/((h^3))*(c^2)).*((E.^2)./(exp(E./(K*Ts))-1));
fun=@(E) (E.*(2*Fs/((h^3)*(c^2))).*((E.^2)./(exp(E./(K*Ts))-1)));
Ps=arrayfun(@(ul)integral(fun,0,ul), E);
figure(1)
plot(E,bs),xlabel('Photon Energy'),ylabel('Spectral Flux Density')
figure(2)
plot(E,Ps),xlabel('Photon Energy'),ylabel('Power Density')
I changed the end value of ‘E’ just to see what the curve did, and discovered that it became asymptotic at about 6.. Change it back if necessary. Also, 250 elements of it are enough to produce a smooth curve.
.
