How do I subtact two or more functions and integrate them?

Hello Everyone

I am trying to obtain a characteristic curve which involves integrating the subtraction of two or more functions . I am facing probleb in the integration part. The below image is the attempted exercise.

I have written a code with my limited knowdge.

clc
clear all

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
q=1.6*10^-16;                 % Electronic Charge
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
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy

E=0:0.0001:4;
V=0:0.0001:4;

bs=((2*Fs)/((h^3))*(c^2)).*((E.^2)./(exp(E./(K*Ts))-1));
b_equilibrium= ((2*Fa)/((h^3))*(c^2)).*((E.^2)./(exp(E./(K*Ta))-1));
b_illumination= ((2*Fa)/((h^3))*(c^2)).*(((E.^2)./(exp((E-(q.*V))./(K*Ta))-1)));


figure(1)
plot(E,bs),xlabel('Photon Energy(eV)'),ylabel('Spectral Flux Density')
figure(2)
plot(E,b_equilibrium),xlabel('Photon Energy(eV)'),ylabel('Dark Flux Density')
figure(3)
plot(E,b_illumination),xlabel('Photon Energy(eV)'),ylabel('Absorbed Flux Density')


fun=@(E) (bs-(b_illumination-b_equilibrium));
Ps=arrayfun(@(ul)integral(fun,Eg,ul), E);



Jv=q.*Ps;
figure(4)
plot(V,Jv),xlabel('Voltage'),ylabel('Current Density')

I am having problem in the integration part. I am having the following error in command window.

Can anyone please help me doing the integration and plot the result. It will be a great help.

Answer 1

[Matt J]

clc
close all
clear all

Fs= 2.16*10^-5*pi;            % Geometrical Factor for sun
Fa=pi;                        % Geometrical Factor for earth
q=1.6*10^-16;                 % Electronic Charge
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
Ta=300;                       % Ambient Temparature
Eg=1.6;                       % BandGap Energy

E=0:0.0001:4;
V=0:0.0001:4;

bs=((2*Fs)/((h^3))*(c^2)).*((E.^2)./(exp(E./(K*Ts))-1));
b_equilibrium= ((2*Fa)/((h^3))*(c^2)).*((E.^2)./(exp(E./(K*Ta))-1));
b_illumination= ((2*Fa)/((h^3))*(c^2)).*(((E.^2)./(exp((E-(q.*V))./(K*Ta))-1)));


figure(1)
plot(E,bs),xlabel('Photon Energy(eV)'),ylabel('Spectral Flux Density')

figure(2)
plot(E,b_equilibrium),xlabel('Photon Energy(eV)'),ylabel('Dark Flux Density')

figure(3)
plot(E,b_illumination),xlabel('Photon Energy(eV)'),ylabel('Absorbed Flux Density')

fun=griddedInterpolant(E,   (bs-(b_illumination-b_equilibrium)));
Ps=arrayfun(@(ul)integral(@(u)fun(u),Eg,ul), E);

Warning: Inf or NaN value encountered.

Jv=q.*Ps;
figure(4)
plot(V,Jv),xlabel('Voltage'),ylabel('Current Density')