I am unable to obtain multiple graphs by varying a parameter in a for loop while solving a Runge-Kutta based problem

While solving a Runge-Kutta problem, I wanted to obtain 3 graphs varying the values of 'phi' in a 'for' loop, but I am unable to obtain it. Please, can someone correct it?

 clear;clf;clc;
%parameters
s=0.01;
t_final=10;
global Re Pr R phi b eps Mn Bi alp rhof cpf kf rhos cps ks h_t
Re= 2.5;
Pr= 6.8;
R= 1;
pp=[0 0.05 0.1];
b=1.0;
eps=0.1;
Mn= 1.0;
Bi= 1.0;
alp= 1;
rhof= 997.1;
cpf= 4179.0;
kf= 0.613;
rhos= 8933.0;
cps= 385.0;
ks= 400;
h_t=2;
phi1= (1-phi)^2.5*((1-phi)+phi*(rhos/rhof));
phi2= (1-phi)+phi*((rhos*cps)/(rhof*cpf));
k    = (ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));
%initial conditions
 t(1)=0;
 h0(1)=1;
 %function handle
 f=@(t,h0) -(1+b)*h0+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*h0^3/3+...
 2*eps*alp*(1-phi)^2.5*(k/(k+Bi*h0))*h0^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
     ((2*k*Bi^2*h0^3*h_t)/(k+Bi*h0)^3+(2*(1+b)*Bi*k*h0^3)/(k+Bi*h0)^3+((1+b)*(3*k+Bi*h0)*k*h0^2)/(k+Bi*h0)^2)*h0^2/2;
 % RK4 loop
 for i=1:numel(pp);
     phi=pp(i);
  for i=1:ceil(t_final/s);
     t(i+1)=t(i)+s;

     k1= f(t(i)       ,  h0(i));
     k2= f(t(i)+0.5*s ,  h0(i)+0.5*k1*s);
     k3= f(t(i)+0.5*s ,  h0(i)+0.5*k2*s);
     k4= f(t(i)+s     ,    h0(i)+k3*s);
     h0(i+1)= h0(i)+s/6*(k1+2*k2+2*k3+k4);
 fprintf('\n%16.3f%16.3f',t,h0);
 plot(t,h0);hold on
  end
  end

Answer 1

[Torsten]

   clear;clf;clc;
  %parameters
  s=0.01;
  t_final=10;
  global Re Pr R phi b eps Mn Bi alp rhof cpf kf rhos cps ks h_t
  Re= 2.5;
  Pr= 6.8;
  R= 1;
  pp=[0 0.05 0.1];
  b=1.0;
  eps=0.1;
  Mn= 1.0;
  Bi= 1.0;
  alp= 1;
  rhof= 997.1;
  cpf= 4179.0;
  kf= 0.613;
  rhos= 8933.0;
  cps= 385.0;
  ks= 400;
  h_t=2;
  for j=1:numel(pp)
    phi=pp(j);
    phi1= (1-phi)^2.5*((1-phi)+phi*(rhos/rhof));
    phi2= (1-phi)+phi*((rhos*cps)/(rhof*cpf));
    k    = (ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));
    %initial conditions
    t(1)=0;
    h(j,1)=1;
    %function handle
    f=@(t,y) -(1+b)*y+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*y^3/3+...
             2*eps*alp*(1-phi)^2.5*(k/(k+Bi*y))*y^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
             ((2*k*Bi^2*y^3*h_t)/(k+Bi*y)^3+(2*(1+b)*Bi*k*y^3)/(k+Bi*y)^3+((1+b)*(3*k+Bi*y)*k*y^2)/(k+Bi*y)^2)*y^2/2;
   % RK4 loop
    for i=1:ceil(t_final/s);
       t(i+1)=t(i)+s;
       k1= f(t(i)       ,  h(j,i));
       k2= f(t(i)+0.5*s ,  h(j,i)+0.5*k1*s);
       k3= f(t(i)+0.5*s ,  h(j,i)+0.5*k2*s);
       k4= f(t(i)+s     ,  h(j,i)+k3*s);
       h(j,i+1)= h(j,i)+s/6*(k1+2*k2+2*k3+k4);
    end
  end
  plot(t,h(1,:),t,h(2,:),t,h(3,:))

And you should preallocate t and h before entering the first for-loop.

Best wishes

Torsten.

[naygarp]

suppose I want the parameter 'h_t' to vary with each iteration of the R-K 4 loop and I used it in this way

    for i=1:ceil(t_final/s);
       t(i+1)=t(i)+s;
       k1= f(t(i)       ,  h(j,i));
       k2= f(t(i)+0.5*s ,  h(j,i)+0.5*k1*s);
       k3= f(t(i)+0.5*s ,  h(j,i)+0.5*k2*s);
       k4= f(t(i)+s     ,  h(j,i)+k3*s);
       h(j,i+1)= h(j,i)+s/6*(k1+2*k2+2*k3+k4);
       h_t=-(1-b)*h(j,i+1);
    end

But I am again getting an error

    Error using  / 
Matrix dimensions must agree.
Error in
>@(t,y)-(1+b)*y+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*y^3/3+2*eps*alp*(1-phi)^2.5*(k/(k+Bi*y))*y^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*((2*k*Bi^2*y^3*h_t)/(k+Bi*y)^3+(2*(1+b)*Bi*k*y^3)/(k+Bi*y)^3+((1+b)*(3*k+Bi*y)*k*y^2)/(k+Bi*y)^2)*y^2/2
(line 38)
       k2= f(t(i)+0.5*s ,
       h(j,i)+0.5*k1*s);

How do I correct it?

[Torsten]

Each time you change a parameter used in the definition of your function, you must update your function handle. Thus after h_t = ..., the line f=@(t,y)... must follow. But are you sure to change h_t after each integration step ?

[naygarp]

yes, 'h_t' would change for each iteration, that is what I'm told.

[naygarp]

will I have to add the function handle 'f=@(t,h)' after in the 2nd for loop

    for i=1:ceil(t_final/s);
       t(i+1)=t(i)+s;
       k1= f(t(i)       ,  h(j,i));
       k2= f(t(i)+0.5*s ,  h(j,i)+0.5*k1*s);
       k3= f(t(i)+0.5*s ,  h(j,i)+0.5*k2*s);
       k4= f(t(i)+s     ,  h(j,i)+k3*s);
       h(j,i+1)= h(j,i)+s/6*(k1+2*k2+2*k3+k4);
       h_t=-(1-b)*h(j,i+1);
       f=@(t,y) -(1+b)*y+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*y^3/3+...
             2*eps*alp*(1-phi)^2.5*(k/(k+Bi*y))*y^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
             ((2*k*Bi^2*y^3*h_t)/(k+Bi*y)^3+(2*(1+b)*Bi*k*y^3)/(k+Bi*y)^3+((1+b)*(3*k+Bi*y)*k*y^2)/(k+Bi*y)^2)*y^2/2;
    end

[Torsten]

It's dirty hack, but I think this will do.

Best wishes

Torsten.