Code does not update the array for each iteration of for loops

My code is not filling updating the temperature functino array (100x100) for each iteration of the for loops. Each column instead is a repeat of the next one. I understand that for the first couple of columns the answer is NaN becausre the solution should be -infinity, but I am not sure why it is simply repeating the previous results. Any help is much appreciated. I haave posted both functions and the script for clarity. Thanks!

function BB = feigR13(num,R,Bi)
x=R-1;
Cc=Bi*x;
p=1.5;
%Preallocation
BB = zeros(1,num);
for j= 1:num
    R11=j*pi/x+0.155*(x)*(j-2)/j^2/(x+2.32)^1.8- 4*(j-1)*(0.062*x/(x+2)^1.56-0.0011)/j^2;
    R12=(j-1/2)*pi/(x)- ((0.305+0.01*x^0.2+0.001*x)/(1+0.5*x)-0.12*(0.2*x^1)/(x^2+1.4)+0.067*x/(x^2+49))/j^1.5;
    lambi=R11*Cc^p/(j*pi+Cc^p)+j*pi*R12/(j*pi+Cc^p);
    lambi=round((lambi*10^10)/(10^10),-20);
    x0 = lambi; 
    dx = pi/x/10;
    for pp=1:6
        UL = (x0 + dx);
        LL = (x0 - dx);
        ML = (UL + LL)/2;
        f0=(-ML*bessely(1,ML*R)+Bi*bessely(0,ML*R))*besselj(0,ML)-(-ML*besselj(1,ML*R)+Bi*besselj(0,ML*R))*bessely(0,ML);
        f1=(-LL*bessely(1,LL*R)+Bi*bessely(0,LL*R))*besselj(0,LL)-(-LL*besselj(1,LL*R)+Bi*besselj(0,LL*R))*bessely(0,LL);
        f2=(-UL*bessely(1,UL*R)+Bi*bessely(0,UL*R))*besselj(0,UL)-(-UL*besselj(1,UL*R)+Bi*besselj(0,UL*R))*bessely(0,UL);
        fp = (f2 - f1)/2/dx;
        fpp = (f1 + f2 - 2*f0)/dx^2;
        h = -f0/fp;
        eps=-fpp*h^2/2/(fp + h*fpp);
        x0 = x0 + h + eps;
        dx = h^2;
        if dx<10^-10
            break
        end
    end
    
    %Assign directly to output variable
    BB(j)=x0;
end
end

clear all
close all
clc
%defining the dimensioned characteristics of the superheater element
A = 10; %systmm accuracy
k = 30; % BTU/hr-ft-F
h = 1000; % BTU/hr-ft^2-F
R1 = 0.137/2; %ft
R2 = 0.167/2; %ft
alpha = 0.0002; %ft^2/s thermal diffusivity for low-carbon steel (AISI 1010)
%function inputs
R = R2/R1;
rv = [0.01:0.01:1];
tv = [0.01:0.01:1];
Bi = (h*R1)/k;
r = length(rv);
t = length(tv);
TR13B01T0 = zeros(r,t);
for ir = 1:r
    for it = 1:t
        TR13B01T0(ir,it) = fDR13B01T0(rv(ir),tv(it),R,Bi,A);
    end
end
plot(TR13B01T0(:,20),tv)
hold on 
plot(TR13B01T0(:,50),tv)
hold on 
plot(TR13B01T0(:,70),tv)
hold on

function [Td,qd]=fDR13B01T0(rv,tv,R,Bi,A)
srv=length(rv);
stv=length(tv);
Temp=zeros(stv,srv); % Preallocating arrays for speed
flux=zeros(stv,srv); % Preallocating arrays for speed
% calculate number of eigenvalues required to obtain solution with accuracy A
mmax1=floor((R-1)/pi*sqrt(A*log(10)/min(tv)));
bet=feigR13(mmax1,R,Bi);% Call the function to get eigenvalues
for ir = 1:srv % begin space loop
    r=rv(ir);
    term1=(Bi*R*log(r))/(1+Bi*R*log(R)); % steady-state temperature solution
    term2=-(Bi*R/r)/(1+Bi*R*log(R)); % steady-state heat flux
    %solution
    for it=1:stv % begin time loop
        t=tv(it);
        mmax=floor((R-1)/pi*sqrt(A*log(10)/t));
        term3 = 0;
        term4 = 0;
        for ii=1:mmax % begin summation loop
            bt=bet(ii);
            V0=-bt*besselj(1,bt*R)+Bi*besselj(0,bt*R);
            S0=-bt*bessely(1,bt*R)+Bi*bessely(0,bt*R);
            Nr = (S0*besselj(0,bt*r)-V0*bessely(0,bt*r));
            L1 = (S0*besselj(1,bt*R)-V0*bessely(1,bt*R));
            Denominator = (Bi^2+bt^2)*besselj(0,bt)*besselj(0,bt)-V0^2;
            Norm = (Denominator)/(besselj(0,bt)*besselj(0,bt));
            term3 = term3 - (pi*pi/2)*R*Bi*(exp(-bt*bt*t))*bt*Nr*L1/Norm;
            Nm = (S0*besselj(1,bt*r)-V0*bessely(1,bt*r));
            term4 = term4 - (pi*pi/2)*Bi*R*exp(-bt*bt*t)*bt*bt*Nm*L1/Norm;
        end
        T(ir) = abs(term1); % Steady-state temperature solution
        Q(ir) = abs(term2); % Steady-state heat flux solution
        Temp(it,ir) = term1 + term3; % Total temperature solution
        flux(it,ir) = term2 + term4; % Total heat flux solution
    end
end
Td=abs(Temp); 
qd=abs(flux);
end

[Walter Roberson]

Re-arranging the code into more usable form:

%defining the dimensioned characteristics of the superheater element
A = 10; %systmm accuracy
k = 30; % BTU/hr-ft-F
h = 1000; % BTU/hr-ft^2-F
R1 = 0.137/2; %ft
R2 = 0.167/2; %ft
alpha = 0.0002; %ft^2/s thermal diffusivity for low-carbon steel (AISI 1010)
%function inputs
R = R2/R1;
rv = [0.01:0.01:1];
tv = [0.01:0.01:1];
Bi = (h*R1)/k;
r = length(rv);
t = length(tv);
TR13B01T0 = zeros(r,t);
for ir = 1:r
    for it = 1:t
        TR13B01T0(ir,it) = fDR13B01T0(rv(ir),tv(it),R,Bi,A);
    end
end
plot(TR13B01T0(:,20),tv)
hold on 
plot(TR13B01T0(:,50),tv)
hold on 
plot(TR13B01T0(:,70),tv)
hold on

function BB = feigR13(num,R,Bi)
    x=R-1;
    Cc=Bi*x;
    p=1.5;
    %Preallocation
    BB = zeros(1,num);
    for j= 1:num
        R11=j*pi/x+0.155*(x)*(j-2)/j^2/(x+2.32)^1.8- 4*(j-1)*(0.062*x/(x+2)^1.56-0.0011)/j^2;
        R12=(j-1/2)*pi/(x)- ((0.305+0.01*x^0.2+0.001*x)/(1+0.5*x)-0.12*(0.2*x^1)/(x^2+1.4)+0.067*x/(x^2+49))/j^1.5;
        lambi=R11*Cc^p/(j*pi+Cc^p)+j*pi*R12/(j*pi+Cc^p);
        lambi=round((lambi*10^10)/(10^10),-20);
        x0 = lambi; 
        dx = pi/x/10;
        for pp=1:6
            UL = (x0 + dx);
            LL = (x0 - dx);
            ML = (UL + LL)/2;
            f0=(-ML*bessely(1,ML*R)+Bi*bessely(0,ML*R))*besselj(0,ML)-(-ML*besselj(1,ML*R)+Bi*besselj(0,ML*R))*bessely(0,ML);
            f1=(-LL*bessely(1,LL*R)+Bi*bessely(0,LL*R))*besselj(0,LL)-(-LL*besselj(1,LL*R)+Bi*besselj(0,LL*R))*bessely(0,LL);
            f2=(-UL*bessely(1,UL*R)+Bi*bessely(0,UL*R))*besselj(0,UL)-(-UL*besselj(1,UL*R)+Bi*besselj(0,UL*R))*bessely(0,UL);
            fp = (f2 - f1)/2/dx;
            fpp = (f1 + f2 - 2*f0)/dx^2;
            h = -f0/fp;
            eps=-fpp*h^2/2/(fp + h*fpp);
            x0 = x0 + h + eps;
            dx = h^2;
            if dx<10^-10
                break
            end
        end
        
        %Assign directly to output variable
        BB(j)=x0;
    end
end

function [Td,qd]=fDR13B01T0(rv,tv,R,Bi,A)
    srv=length(rv);
    stv=length(tv);
    Temp=zeros(stv,srv); % Preallocating arrays for speed
    flux=zeros(stv,srv); % Preallocating arrays for speed
    % calculate number of eigenvalues required to obtain solution with accuracy A
    mmax1=floor((R-1)/pi*sqrt(A*log(10)/min(tv)));
    bet=feigR13(mmax1,R,Bi);% Call the function to get eigenvalues
    for ir = 1:srv % begin space loop
        r=rv(ir);
        term1=(Bi*R*log(r))/(1+Bi*R*log(R)); % steady-state temperature solution
        term2=-(Bi*R/r)/(1+Bi*R*log(R)); % steady-state heat flux
        %solution
        for it=1:stv % begin time loop
            t=tv(it);
            mmax=floor((R-1)/pi*sqrt(A*log(10)/t));
            term3 = 0;
            term4 = 0;
            for ii=1:mmax % begin summation loop
                bt=bet(ii);
                V0=-bt*besselj(1,bt*R)+Bi*besselj(0,bt*R);
                S0=-bt*bessely(1,bt*R)+Bi*bessely(0,bt*R);
                Nr = (S0*besselj(0,bt*r)-V0*bessely(0,bt*r));
                L1 = (S0*besselj(1,bt*R)-V0*bessely(1,bt*R));
                Denominator = (Bi^2+bt^2)*besselj(0,bt)*besselj(0,bt)-V0^2;
                Norm = (Denominator)/(besselj(0,bt)*besselj(0,bt));
                term3 = term3 - (pi*pi/2)*R*Bi*(exp(-bt*bt*t))*bt*Nr*L1/Norm;
                Nm = (S0*besselj(1,bt*r)-V0*bessely(1,bt*r));
                term4 = term4 - (pi*pi/2)*Bi*R*exp(-bt*bt*t)*bt*bt*Nm*L1/Norm;
            end
            T(ir) = abs(term1); % Steady-state temperature solution
            Q(ir) = abs(term2); % Steady-state heat flux solution
            Temp(it,ir) = term1 + term3; % Total temperature solution
            flux(it,ir) = term2 + term4; % Total heat flux solution
        end
    end
    Td=abs(Temp); 
    qd=abs(flux);
end

[Walter Roberson]

Re-arranging the plotting as you had reversed x and y coordinates (but it does not make any difference to the overall problem)

%defining the dimensioned characteristics of the superheater element
A = 10; %systmm accuracy
k = 30; % BTU/hr-ft-F
h = 1000; % BTU/hr-ft^2-F
R1 = 0.137/2; %ft
R2 = 0.167/2; %ft
alpha = 0.0002; %ft^2/s thermal diffusivity for low-carbon steel (AISI 1010)
%function inputs
R = R2/R1;
rv = [0.01:0.01:1];
tv = [0.01:0.01:1];
Bi = (h*R1)/k;
r = length(rv);
t = length(tv);
TR13B01T0 = zeros(r,t);
for ir = 1:r
    for it = 1:t
        TR13B01T0(ir,it) = fDR13B01T0(rv(ir),tv(it),R,Bi,A);
    end
end
plot(tv, TR13B01T0(:,20), 'displayname', string(tv(20)))
hold on 
plot(tv, TR13B01T0(:,50), 'displayname', string(tv(50)))
hold on 
plot(tv, TR13B01T0(:,70), 'displayname', string(tv(70)))
hold on
legend show

max((TR13B01T0(:,20) - TR13B01T0(:,50)).')

ans = 0

function BB = feigR13(num,R,Bi)
    x=R-1;
    Cc=Bi*x;
    p=1.5;
    %Preallocation
    BB = zeros(1,num);
    for j= 1:num
        R11=j*pi/x+0.155*(x)*(j-2)/j^2/(x+2.32)^1.8- 4*(j-1)*(0.062*x/(x+2)^1.56-0.0011)/j^2;
        R12=(j-1/2)*pi/(x)- ((0.305+0.01*x^0.2+0.001*x)/(1+0.5*x)-0.12*(0.2*x^1)/(x^2+1.4)+0.067*x/(x^2+49))/j^1.5;
        lambi=R11*Cc^p/(j*pi+Cc^p)+j*pi*R12/(j*pi+Cc^p);
        lambi=round((lambi*10^10)/(10^10),-20);
        x0 = lambi; 
        dx = pi/x/10;
        for pp=1:6
            UL = (x0 + dx);
            LL = (x0 - dx);
            ML = (UL + LL)/2;
            f0=(-ML*bessely(1,ML*R)+Bi*bessely(0,ML*R))*besselj(0,ML)-(-ML*besselj(1,ML*R)+Bi*besselj(0,ML*R))*bessely(0,ML);
            f1=(-LL*bessely(1,LL*R)+Bi*bessely(0,LL*R))*besselj(0,LL)-(-LL*besselj(1,LL*R)+Bi*besselj(0,LL*R))*bessely(0,LL);
            f2=(-UL*bessely(1,UL*R)+Bi*bessely(0,UL*R))*besselj(0,UL)-(-UL*besselj(1,UL*R)+Bi*besselj(0,UL*R))*bessely(0,UL);
            fp = (f2 - f1)/2/dx;
            fpp = (f1 + f2 - 2*f0)/dx^2;
            h = -f0/fp;
            eps=-fpp*h^2/2/(fp + h*fpp);
            x0 = x0 + h + eps;
            dx = h^2;
            if dx<10^-10
                break
            end
        end
        
        %Assign directly to output variable
        BB(j)=x0;
    end
end

function [Td,qd]=fDR13B01T0(rv,tv,R,Bi,A)
    srv=length(rv);
    stv=length(tv);
    Temp=zeros(stv,srv); % Preallocating arrays for speed
    flux=zeros(stv,srv); % Preallocating arrays for speed
    % calculate number of eigenvalues required to obtain solution with accuracy A
    mmax1=floor((R-1)/pi*sqrt(A*log(10)/min(tv)));
    bet=feigR13(mmax1,R,Bi);% Call the function to get eigenvalues
    for ir = 1:srv % begin space loop
        r=rv(ir);
        term1=(Bi*R*log(r))/(1+Bi*R*log(R)); % steady-state temperature solution
        term2=-(Bi*R/r)/(1+Bi*R*log(R)); % steady-state heat flux
        %solution
        for it=1:stv % begin time loop
            t=tv(it);
            mmax=floor((R-1)/pi*sqrt(A*log(10)/t));
            term3 = 0;
            term4 = 0;
            for ii=1:mmax % begin summation loop
                bt=bet(ii);
                V0=-bt*besselj(1,bt*R)+Bi*besselj(0,bt*R);
                S0=-bt*bessely(1,bt*R)+Bi*bessely(0,bt*R);
                Nr = (S0*besselj(0,bt*r)-V0*bessely(0,bt*r));
                L1 = (S0*besselj(1,bt*R)-V0*bessely(1,bt*R));
                Denominator = (Bi^2+bt^2)*besselj(0,bt)*besselj(0,bt)-V0^2;
                Norm = (Denominator)/(besselj(0,bt)*besselj(0,bt));
                term3 = term3 - (pi*pi/2)*R*Bi*(exp(-bt*bt*t))*bt*Nr*L1/Norm;
                Nm = (S0*besselj(1,bt*r)-V0*bessely(1,bt*r));
                term4 = term4 - (pi*pi/2)*Bi*R*exp(-bt*bt*t)*bt*bt*Nm*L1/Norm;
            end
            T(ir) = abs(term1); % Steady-state temperature solution
            Q(ir) = abs(term2); % Steady-state heat flux solution
            Temp(it,ir) = term1 + term3; % Total temperature solution
            flux(it,ir) = term2 + term4; % Total heat flux solution
        end
    end
    Td=abs(Temp); 
    qd=abs(flux);
end

Answer 1

[Walter Roberson]

        lambi=round((lambi*10^10)/(10^10),-20);

Notice the round() covers all of (lambi*10^10)/(10^10) as its first parameter. But algebraically the multiplication and division cancel out, which would lead to round(lambi,-20) which is a request to round lambi to the nearest 10^20 . The existing code only differs from that to the extent that if lambi is sufficiently large then lambi*10^10 might not be exactly representable, so potentially the (lambi*10^10)/(10^10) changes the round-off of lambi.. or that lambi*10^10 might overflow to infinity. An example of a value that changes in going through *10^10/10^10 but is small enough to be affected by the round(-20), is lambi = 6.1479804705177653e+69

At the very least, that section needs to be heavily commented to describe why it is important to do that multiplying and diviging by 10^10 and that rounding to the nearest 10^20 .

The actual values that go into calculating lambi are relative constants -- the same for every pairs of r and t values. So before the odd rounding, lambi will always be roughly 7 2/3. There are values near that, such as 7.2664788935685998 that come out slightly different with the *10^10/10^10 process... but all values in that range round to 0 to the nearest 10^20 .

So... your lambi is rounding to 0. And then when you step through the UL, LL, ML calculation, ML will come out as 0. You have bessely(0,ML) but that is -inf ... so you have calculations involving infinities, and those will typically generate NaN.

The reason why you do not always generate NaN for all of your iterations, is that current value of t makes a difference in calculationg mmax1, and mmax1 is used as the upper limit of your for loop. But with t values beyond roughly 0.11 then mmax1 comes out as 0, so your loop is 1:0 which causes it not to be done at all, leaving the arrays still at the zero() that you initialized them to.

Basically your code is calculating whether t is more than some number greater than roughly 0.11, and if so then return 0, and otherwise return NaN.

[Matthew Palermo]

Gotcha, I was overlooking that. I got this code from a scientific paper, but I have had issues using the authors' codes before.

Thanks!