I an not getting multiple graphs after running two 'for' loops for different parameters solving ode using 'bvp4c'

The equations that I have to solve are :

f"=g(g^2+gamma^2)/(g^2+lambda*gamma^2) ------ (1)

g'= 1/3f'^2-2/3(f*f") ------------------------(2)

t"+4*Rd*t"/3+ 2*Pr*f*t'/3+ Nb*t'*p'+Nt*(t')^2= 0------(3)

p"+(2*Lew*f*p')/3+ Nt*t"/Nb= 0 ------------------------(4)

Where 'lambda', 'gamma', 'Rd', 'Pr', 'Lew', 'Nb', 'Nt' are some parameters

where the boundary conditions are : f=0, f'= delta, t=1, p=1 at y=0 f'-->0, t-->0, p-->0 as y----> infinity

I tried to solve these equations using 'bvp4c'

I am not getting multiple graphs even after running two 'for' loops for different values of parameters 'lambda', 'Rd'

Please help.

function sol= proj
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt delta
delta=1;
gama=1;
Pr=2;
Lew=2; 
Nb=1;
Nt=2;
pp=[0.5:0.5:1.5];
qq=[0 1 2];
figure(1)
plot(2,1);hold on
solinit= bvpinit([0:0.01:10],[1,0,0,0,0,0,0]);
for i= 1:numel(pp)
    k=1;
    for lambda=pp(i)
    Rd=qq(k);
sol= bvp4c(@projfun,@projbc,solinit);
k=k+1;
solinit= sol;
    end  
end
plot(2,1);plot(sol.x,sol.y(4,:));hold on 
end
function f= projfun(x,y)
global lambda gama Pr Rd Lew Nb Nt
f= [y(2)
    y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2)
    y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))
    y(5)
    -((2*Pr*y(1)*y(5))/3 + Nb*y(5)*y(7) + Nt*y(5)^2)/(1+4*Rd/3)
    y(7)
    -(2*Lew*y(1)*y(7))+ Nt*((2*Pr*y(1)*y(5))/3 + Nb*y(5)*y(7) + Nt*y(5)^2)/Nb*(1+4*Rd/3)];
end
function res= projbc(ya,yb)
global delta
res= [ya(1); ya(2)-delta; ya(4)-1.0; ya(6)-1.0; yb(2); yb(4); yb(6)];
end

[naygarp]

Shouldn't I be getting 9 graphs in one figure here?

[MINATI]

Dear naygarp have u drawn the figs u asked from W. A. Khan & I. Pop paper for Nt and Nb perhaps. If yes can u pl. share me the matlab program at [e-mail] Thanks in advance

[naygarp]

I could only extract the values of Reduced Nusslet's no(-theta'(eta)). and Reduced Sherwood's no.(-phi'(eta)) and plot it manually, extracting and putting those values vs the values of Nb or Nt. It takes a while to copy and paste all those values.

Answer 1

[Walter Roberson]

Your code does not ask for multiple graphs.

You start with

    plot(2,1);hold on

which designates a single point. You do not specify at least two points and you do not specify a marker, so nothing shows up on the display when you do this. Mostly it has an effect on the automatic axes boundary selection, otherwise no effect on the display.

You then have

    for i= 1:numel(pp)
        k=1;
        for lambda=pp(i)
            Rd=qq(k);
            sol= bvp4c(@projfun,@projbc,solinit);
            k=k+1;
            solinit= sol;
        end
    end

which kind of looks like a double nested for loop, but the "for lambda=pp(i)" is just a fancy way of saying "lambda=pp(i)" there since i is not a vector. So anyhow, you have a for i loop there, and there are no graphics within the loop. I am not clear as to why you are feeding solinit back through several rounds, but I guess there might be a reason for that.

After the loop you have

    plot(2,1);plot(sol.x,sol.y(4,:));hold on

the plot(2,1) plots on top of the previous plot(2,1) -- except, as noted before, since that is only a single point and you did not specify a marker, that doesn't change the display (except perhaps the xlim and ylim chosen.)

The plot(sol.x,sol.y(4,:)) part selects the 4th output variable out of the 7 output variables and plots a line with just that one variable. So you get one line out of that.

... and that's the end of your program. You plotted the same invisible point twice, and you plot one line. The best you could hope for with that is the line together with a marker at (2,1) -- at most 2 outputs, not 9.

[naygarp]

Thanks for your reply. What I want are graphs with a variation of say any two parameters of my equation like the picture of the graph I have attached herewith. For variation of one parameter a simple for loop does the job. Now for variation of two parameters shouldn't I have to use two for loops.

Please if possible could you show what modifications do I have to make.

For variation of one parameter 'lambda' I have used the similar code code and I get the graphs

function sol= proj
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt delta
delta=1;
gama=1;
Pr=2;
Rd=1;
Lew=2; 
Nb=1;
Nt=2;
pp=[0.5:0.5:1.5];
figure(1)
plot(2,1);hold on
solinit= bvpinit([0:0.01:10],[1,0,0,0,0,0,0]);
for lambda= pp
sol= bvp4c(@projfun,@projbc,solinit);
solinit= sol;
plot(2,1);plot(sol.x,sol.y(4,:));hold on
end
end
function f= projfun(x,y)
global lambda gama Pr Rd Lew Nb Nt
f= [y(2)
  y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2)
  y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))
  y(5)
  -((2*Pr*y(1)*y(5))/3 + Nb*y(5)*y(7) + Nt*y(5)^2)/(1+4*Rd/3)
  y(7)
  -(2*Lew*y(1)*y(7))+ Nt*((2*Pr*y(1)*y(5))/3 + Nb*y(5)*y(7) + Nt*y(5)^2)/Nb*(1+4*Rd/3)];
end
function res= projbc(ya,yb)
global delta
res= [ya(1); ya(2)-delta; ya(4)-1.0; ya(6)-1.0; yb(2); yb(4); yb(6)];
end

[naygarp]

How can I show which curve is for which parameter value in the plot box. Is the 'legend' command used to do that?

[Walter Roberson]

Yes, use legend() for that.

[naygarp]

If I use like the following, I get which curve is for which parameter, Is it showing accurately or not please could you check it?

    plot(2,1);plot(sol.x,sol.y(4,:));hold on  
legend('1st value','2nd value','3rd value');

[Walter Roberson]

plot() of a single point with no marker is almost always wrong (except that plot(nan,nan) or plot(inf,inf) can have uses in customizing legends.)

The hold on goes after the first plot(), not after the second plot()