How can I solve the function numerically for changing Prandtl Numbers ? I couldn't handle this with a for loop inside of a function.

First of all, Hi to everyone. I hope you're okay in these days.

I would like to have a temperature distrubution for 5 different Prandtl numbers but I failed. The for loop that I created is not working. The function only takes the last index, so the last prandtl number is taken, and calculations have done according to this Prandtl number.

But I want to have the solution for different Prandtl Numbers and plot all of the solutions according to eta values.

My question is, the most important one here is;

1.How can I have solutions for different Prandtl numbers ?+++

and less importantly ,how can I plot all the solutions in a single graph ? ---

Have a nice and healty days.

This is my code:

function [Tstar]=Tstar_bools_(a,b,N)
if nargin<3,a=0;b=5; N=1600; end
eta = (0:0.2:b);
numberofeta=length(eta);
A = zeros(numberofeta,2);
B=ones(numberofeta,2);

for k = 1:numberofeta
    
    h=(eta(k)-a)/N;
    sum=(h/90)*f(a)+(h/90)*f(eta(k));
    for i=1:N-1
        x(i)=a+(i*h);
        sum=sum+((2*h)/180)*(7*f(x(i))+32*f(x(i)+h/4)+12*f(x(i)+(2*h)/4)+32*f(x(i)+(3*h)/4)+7*f(x(i)+h));        
        A(k,1)=eta(k);
        A(k,2)=sum; 
        B(1:k,2)=A(numberofeta,2);    
    end
    
end
A
% B

Tstar = A./B;
% Tstar;
% plot(Tstar(:,1), Tstar(:,2:end),'ro-')
end

function fx = f(x)
    Pr=[0.7 1 3 7 13];
    for i = 1:1:5
    fx = (2*10^-5*x - 0.0011*x + 0.0163*x - 0.0843*x + 0.0824*x + 0.3177).^(Pr(i));
    end
end

% function fx = f(x)
%     fx = (2*10^-5*x - 0.0011*x + 0.0163*x - 0.0843*x + 0.0824*x + 0.3177)^(13);
% end

Answer 1

[Alan Stevens]

Like this?

a = 0; b = 5; N = 1600;
eta = (0:0.2:b);
numberofeta=length(eta);
A = zeros(numberofeta,6); %%%%%%%%
B=ones(numberofeta,6);    %%%%%%%%

for k = 1:numberofeta
    
    h=(eta(k)-a)/N;
    sum=(h/90)*f(a)+(h/90)*f(eta(k));
    for i=1:N-1
        x(i)=a+(i*h);
        sum=sum+((2*h)/180)*(7*f(x(i))+32*f(x(i)+h/4)+12*f(x(i)+(2*h)/4)+32*f(x(i)+(3*h)/4)+7*f(x(i)+h));        
        A(k,1)=eta(k);
        A(k,2:6)=sum; %%%%%%%%%%
        B = A(numberofeta,:); %%%%%%% 
    end
    
end

Tstar = A./B;
Tstar;
plot(Tstar(:,1), Tstar(:,2:end),'o-')

function fx = f(x)
    Pr=[0.7 1 3 7 13];
    for i = 1:1:5
    fx(i) = (2*10^-5*x - 0.0011*x + 0.0163*x - 0.0843*x + 0.0824*x + 0.3177).^(Pr(i)); %%%%%%
    end
end

[Baris Ecevit]

Yes exactly like this. I'm very happy. I really thank you for this. I check the code and everything is fine really but in the Matrix B, the first column has to be remain "1" in order not to change the first column of A matrix .because the first column is the eta values and they're the x values of the plot. So, The x axis should be up to 5. I'm trying to fix this now. I mean in this line of code ;

B = A(numberofeta,:); %%%%%%%

B matrix has to change except the column number 1. Ijust want to inform you about this. You did great work for me and I appreciate it :D

Have a nice and healty days.

Edit to my last comment:

I solved the problem. I changed the first column by "1" by adding a line of code.