I want 4 solid lines for Fr paramter . but im getting single line why?

function S_a
clc
clear All
                

                        %%%%note%%%%%%%%%


%%%%%%%%%%%%%%%%% decreasing behaviour with Fr=0.4,0.8,1.2,1.6,2  %%%%%%%%%%%



global  Fr Pr f2 f1 F2 shi2 shi3 lambda zeta  Sk u Re Kf B  A H1 H2 H3 C1 C2 C3 P1 P2 P3 K1 K2 K3  Kb n Z  FR


FR = 0.4:0.4:2
for i=1:numel(FR) 
Fr = FR(i); %variation paramter
end
Pr = 6.2 ; 
B = 0.9;  
F2 = 0:0.01:0.15
for i=1:numel(F2) 
f2 = F2(i);
    f1 = 0.05 ; 
    lambda = 0.4; 
    zeta=0.5; 
    A= 0.7;  
    n=3.7;  
    Re=0.3; 
    
    
    C1=765;
    P1=3970;
    K1=40;
    %%%%%%%%%%%%%%%%%%%%%%%
    C2=385; % specific heat
    P2=8933; % density
    K2=400; % thermal conductivity
    
    %%%%%%%%%%% %%%%%%%%%%%%
    C3=4180; % specific heat
    P3=997.1; % density
    K3=0.6071; % thermal conductivity
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%
    Z=zeta+B;
    H1=P1*C1;       % pho*cp nanoparticle 1
    H2=P2*C2;       % pho*cp nanoparticle 2
    H3=P3*C3;       % pho*cp base fluid
    
    Kb = (K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)); % khnf/kbf
    
    Kf=((K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)))*((K1 + (n - 1)*K3 - f1*(n - 1)*(K3 - K1))/(K1 + (n - 1)*K3 + f1*(K3 - K1)));
    % khnf/kf
    
    shi2=(((1 - f1)^2.5)*((1 - f2)^2.5))*(((1 - f1)*(1 - f2))+f1*(P1/P3)) + f2*(P2/P3)  ;  
    % shi2 in  velocity equation
    
    
    shi3=(1 - f2)*((1 - f1) + f1*((H1)/(H3))) + f2*((H2)/(H3)); 
    
    u=1/((1-f1)^2.5);
    
    
    options = bvpset('RelTol',1e-6,'Stats','off') ;  
    solinit = bvpinit(linspace(0,5,10),[1; 0 ;0 ;0 ;0; 0]);
    sol = bvp4c(@odes,@bcs,solinit,options);
    Sk(i)=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
    
end
plot(F2,Sk,'-y')
%--------------------------------------------------------------------------
    function dydx = odes(~,y)
        % global zeta Pr psi2 psi3 k B  Fr A
        dydx = zeros(6,1);
        dydx(1) = y(2); % y(1)=f
        dydx(2) = y(3);
        dydx(3) = y(4);
        dydx(4) = -(2/(zeta+B)*y(4))+(1/((zeta+B)^2)*y(3))-(1/((zeta+B)^3)*y(2))+2*lambda*(y(3)+(1/(zeta+B)*y(2)))+shi2*(-((B/(zeta+B)^2)*y(1)*y(3))-((B/(zeta+B)*y(1)*y(4)))+(B/((zeta+B)^3)*y(1)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(3))+2*Fr*(2*y(2)*y(3)+((1/(zeta+B))*y(2)*y(2))));
        dydx(5) = y(6); % y(5)=theta
        dydx(6) = -(1/(zeta+B))*y(6)-((Pr*shi3*B)/(Kb*(zeta+B)))*(y(1)*y(6)-A*y(2)*y(5));% k=khnf/kbf
        %--------------------------------------------------------------------------
    end
    function res = bcs(ya,yb)
        res = [ya(1); ya(2)-1;ya(5)-1; yb(2); yb(3); yb(5)];
    end

%%%%%%%%%%%%2nd%%%%%%%%%%%%%

Fr = 0.8;  
Pr = 6.2 ; 
B = 0.9;  
F2 = 0:0.01:0.15
for i=1:numel(F2) 
f2 = F2(i);
    f1 = 0.05 ; 
    lambda = 0.4; 
    zeta=0.5; 
    A= 0.7;  
    n=3.7;  
    Re=0.3; 
    
    
    C1=765;
    P1=3970;
    K1=40;
    %%%%%%%%%%%%%%%%%%%%%%%
    C2=385; % specific heat
    P2=8933; % density
    K2=400; % thermal conductivity
    
    %%%%%%%%%%% %%%%%%%%%%%%
    C3=4180; % specific heat
    P3=997.1; % density
    K3=0.6071; % thermal conductivity
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%
    Z=zeta+B;
    H1=P1*C1;       % pho*cp nanoparticle 1
    H2=P2*C2;       % pho*cp nanoparticle 2
    H3=P3*C3;       % pho*cp base fluid
    
    Kb = (K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)); % khnf/kbf
    
    Kf=((K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)))*((K1 + (n - 1)*K3 - f1*(n - 1)*(K3 - K1))/(K1 + (n - 1)*K3 + f1*(K3 - K1)));
    % khnf/kf
    
    shi2=(((1 - f1)^2.5)*((1 - f2)^2.5))*(((1 - f1)*(1 - f2))+f1*(P1/P3)) + f2*(P2/P3)  ;  
    % shi2 in  velocity equation
    
    
    shi3=(1 - f2)*((1 - f1) + f1*((H1)/(H3))) + f2*((H2)/(H3)); 
    
    u=1/((1-f1)^2.5);
    
    
    options1 = bvpset('RelTol',1e-6,'Stats','off') ;  
    solinit1 = bvpinit(linspace(0,5,10),[1; 0 ;0 ;0 ;0; 0]);
    sol = bvp4c(@odes1,@bcs1,solinit1,options1);
    Sk(i)=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
    
end
plot(F2,Sk,'--r')
%--------------------------------------------------------------------------
    function dydx = odes1(~,y)
        % global zeta Pr psi2 psi3 k B  Fr A
        dydx = zeros(6,1);
        dydx(1) = y(2); % y(1)=f
        dydx(2) = y(3);
        dydx(3) = y(4);
        dydx(4) = -(2/(zeta+B)*y(4))+(1/((zeta+B)^2)*y(3))-(1/((zeta+B)^3)*y(2))+2*lambda*(y(3)+(1/(zeta+B)*y(2)))+shi2*(-((B/(zeta+B)^2)*y(1)*y(3))-((B/(zeta+B)*y(1)*y(4)))+(B/((zeta+B)^3)*y(1)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(3))+2*FR(i)*(2*y(2)*y(3)+((1/(zeta+B))*y(2)*y(2))));
        dydx(5) = y(6); % y(5)=theta
        dydx(6) = -(1/(zeta+B))*y(6)-((Pr*shi3*B)/(Kb*(zeta+B)))*(y(1)*y(6)-A*y(2)*y(5));% k=khnf/kbf
        %--------------------------------------------------------------------------
    end
  
    function res = bcs1(ya,yb)
        res = [ya(1); ya(2)-1;ya(5)-1; yb(2); yb(3); yb(5)];
    end

end

[Torsten]

Why 4 lines ? FR has 5 elements.

[farooq omer]

my bad,for 5 values of Fr . thank you for your response

Answer 1

[Torsten]

S_a()

function S_a               
                          %%%%note%%%%%%%%%
  %%%%%%%%%%%%%%%%% decreasing behaviour with Fr=0.4,0.8,1.2,1.6,2  %%%%%%%%%%%

  FR = 0.4:0.4:2;
  F2 = 0:0.01:0.15;
  Pr = 6.2 ; 
  B = 0.9;
  f1 = 0.05 ; 
  lambda = 0.4; 
  zeta=0.5; 
  A= 0.7;  
  n=3.7;  
  Re=0.3; 
    
  C1=765;
  P1=3970;
  K1=40;
  %%%%%%%%%%%%%%%%%%%%%%%
  C2=385; % specific heat
  P2=8933; % density
  K2=400; % thermal conductivity
    
  %%%%%%%%%%% %%%%%%%%%%%%
  C3=4180; % specific heat
  P3=997.1; % density
  K3=0.6071; % thermal conductivity
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%
  Z=zeta+B;
  H1=P1*C1;       % pho*cp nanoparticle 1
  H2=P2*C2;       % pho*cp nanoparticle 2
  H3=P3*C3;       % pho*cp base fluid

  hold on

  for j=1:numel(FR) 
    Sk = zeros(size(F2)); 
    Fr = FR(j); %variation paramter
    for i=1:numel(F2) 
      f2 = F2(i);
    
      Kb = (K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)); % khnf/kbf
    
      Kf=((K2 + (n - 1)*K3 - f2*(n - 1)*(K3 - K2))/(K2 + (n - 1)*K3 + f2*(K3 - K2)))*((K1 + (n - 1)*K3 - f1*(n - 1)*(K3 - K1))/(K1 + (n - 1)*K3 + f1*(K3 - K1)));
      % khnf/kf
    
      shi2=(((1 - f1)^2.5)*((1 - f2)^2.5))*(((1 - f1)*(1 - f2))+f1*(P1/P3)) + f2*(P2/P3)  ;  
      % shi2 in  velocity equation
    
    
      shi3=(1 - f2)*((1 - f1) + f1*((H1)/(H3))) + f2*((H2)/(H3)); 
    
      u=1/((1-f1)^2.5);
    
    
      options = bvpset('RelTol',1e-6,'Stats','off') ;  
      solinit = bvpinit(linspace(0,5,10),[1; 0 ;0 ;0 ;0; 0]);
      sol = bvp4c(@odes,@bcs,solinit,options);
      Sk(i)=u*1/((Re)^0.5)*(sol.y(3,1)-sol.y(2,1)/B);
    
    end
    plot(F2,Sk,'-y')
  end
  hold off
  %--------------------------------------------------------------------------
    function dydx = odes(~,y)
        % global zeta Pr psi2 psi3 k B  Fr A
        dydx = zeros(6,1);
        dydx(1) = y(2); % y(1)=f
        dydx(2) = y(3);
        dydx(3) = y(4);
        dydx(4) = -(2/(zeta+B)*y(4))+(1/((zeta+B)^2)*y(3))-(1/((zeta+B)^3)*y(2))+2*lambda*(y(3)+(1/(zeta+B)*y(2)))+shi2*(-((B/(zeta+B)^2)*y(1)*y(3))-((B/(zeta+B)*y(1)*y(4)))+(B/((zeta+B)^3)*y(1)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(2))+(3*B/((zeta+B)^2)*y(2)*y(3))+2*Fr*(2*y(2)*y(3)+((1/(zeta+B))*y(2)*y(2))));
        dydx(5) = y(6); % y(5)=theta
        dydx(6) = -(1/(zeta+B))*y(6)-((Pr*shi3*B)/(Kb*(zeta+B)))*(y(1)*y(6)-A*y(2)*y(5));% k=khnf/kbf
        %--------------------------------------------------------------------------
    end
    function res = bcs(ya,yb)
        res = [ya(1); ya(2)-1;ya(5)-1; yb(2); yb(3); yb(5)];
    end
end