% Eq. E.19 simulation

t_a=0.1;
t=[0:0.001:t_a]; %t_az=0.1;

y= pi*exp(2.*t);


% We compute max value of second derivative of E_0(-t) by taking maximum of Positive terms terms_pos_max
% and Minimum of negative terms terms_neg_min.

y_min=min(y); y_max=max(y);

sum_n_1_9=0;

for n=1:2:9, %100,

term_1_pos_max = exp(-n.^2.*y_min/4).*y_max.^(1/4).*(n.^4/16);
term_1_neg_min = exp(-n.^2.*y_max/4).*y_min.^(1/4).*(-n.^2./(8.*y_max) -3./(16.*y_max.^2) );

term_2_pos_max = - exp(-n.^2.*y_min).*y_max.^(1/4).*( -n.^2./(2.*y_min) -3./(16.*y_min.^2));
term_2_neg_min = - exp(-n.^2.*y_max).*y_min.^(1/4).*(n.^4 );

term_3_pos_max = - exp(-(n+1).^2.*y_min/4).*y_max.^(1/4).*(-(n+1).^2./(8.*y_min) -3./(16.*y_min.^2));
term_3_neg_min = - exp(-(n+1).^2.*y_max/4).*y_min.^(1/4).*((n+1).^4/16 );

term_4_pos_max = exp(-(n+1).^2.*y_min).*y_max.^(1/4).*((n+1).^4);
term_4_neg_min = exp(-(n+1).^2.*y_max).*y_min.^(1/4).*(-(n+1).^2./(2.*y_max) -3./(16.*y_max.^2) );

terms_pos_max = term_1_pos_max + term_2_pos_max + term_3_pos_max + term_4_pos_max; 
terms_neg_min = term_1_neg_min + term_2_neg_min + term_3_neg_min + term_4_neg_min; 

%terms_pos_max*(pi).^(1/4)
%terms_neg_min*(pi).^(1/4)

n
(terms_pos_max+terms_neg_min)*(pi).^(-1/4)

sum_n_1_9 = sum_n_1_9 + (terms_pos_max+terms_neg_min)*(pi).^(-1/4);

end

disp('Sum for n=1,3,5,7,9.  .....')
sum_n_1_9 % -0.009792940445532

%sqrt(pi)/2(sqrt(K) erfc(11 sqrt(K)) for n=11,13,...

y_a = pi*exp(2.*t_a); %3.8371

K= y_min/4 - 0.1;
term_1 = (200/16)*  (sqrt(pi)/(2*sqrt(K))) *  erfc(11*sqrt(K));

K= y_min;
term_2 = (3/(16*y_min.^2)) *  (sqrt(pi)/(2*sqrt(K))) *  erfc(11*sqrt(K));

K= y_min - 0.1;
term_3 = (10/(2*y_min)) *  (sqrt(pi)/(2*sqrt(K))) *  erfc(11*sqrt(K));

K= y_min/4;
term_4 = (3/(16*y_min.^2)) *  (sqrt(pi)/(2*sqrt(K))) *  erfc(11*sqrt(K));

K= y_min/4 - 0.1;
term_5 = ( 1.1.^2*10/(8*y_min)) *  ( sqrt(pi)/(2*sqrt(K))) *  erfc(11*sqrt(K));

K= y_min - 0.1;
term_6 = 200*(1.1).^4 *  (sqrt(pi)/(2*sqrt(K))) *  erfc(11*sqrt(K));

sum_n_11_inf = (pi).^(-1/4) * y_max.^(1/4) * (term_1 + term_2 + term_3 + term_4 + term_5 + term_6) ;

disp('Sum for n=11 to infinity....')
sum_n_11_inf % 8.643121748152311e-37




