For such series, if the sum is convergent, it is useful to approximate it with finite terms instead of using infinite series. You can use sum() and integral() functions to implement this equation. See the following code to see how it can be done
xf = 1;
xe = 2;
xd = 3;
sum1_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_);
sum2_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_).*sin(n.*pi.*xf./xe).*cos(n.*pi.*xd.*xf./xe)./(n.*pi.*xf./xe);
sum1 = @(tda_,N) sum(sum1_terms(tda_,1:N));
sum2 = @(tda_,N) sum(sum2_terms(tda_,1:N));
dPwd = @(tda_,N) (1+2*sum1(tda_,N)).*(1+2*sum2(tda_,N));
Pwd = @(tda, N) integral(@(tda_) dPwd(tda_, N), 0, tda, 'ArrayValued', 1);
N_value = 100; % 100 terms will be used in calculations
tda_value = linspace(0,0.01,100);
Pwd_values = zeros(size(tda_value));
for i=1:numel(tda_value)
Pwd_values(i) = Pwd(tda_value(i), N_value);
end
plot(tda_value, Pwd_values)
Following code is a faster version of the above code (about 3x faster), but requires a bit deeper understanding of MATLAB functions
xf = 1;
xe = 2;
xd = 3;
sum1_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_);
sum2_terms = @(tda_,n) exp(-4.*n.^2.*pi^2.*tda_).*sin(n.*pi.*xf./xe).*cos(n.*pi.*xd.*xf./xe)./(n.*pi.*xf./xe);
sum1 = @(tda_,N) sum(sum1_terms(tda_,(1:N).'));
sum2 = @(tda_,N) sum(sum2_terms(tda_,(1:N).'));
dPwd = @(tda_,N) (1+2*sum1(tda_,N)).*(1+2*sum2(tda_,N));
Pwd = @(tda, N) integral(@(tda_) dPwd(tda_, N), 0, tda);
N_value = 100; % 100 terms will be used in calculations
tda_value = linspace(0,0.01,100);
Pwd_values = zeros(size(tda_value));
for i=1:numel(tda_value)
Pwd_values(i) = Pwd(tda_value(i), N_value);
end
plot(tda_value, Pwd_values)
