plotting a graph of multivariable function

I formulated a matlab code for obtaining a curve and it is giving results but inside a loop for y1 and y2. The code is taking only last values, i.e., y1 and y2=1. however, I want to include the values of y1 and y2 from 0 to 1. How can I include it.

The code is written as:

clear all
clc
format longEng
syms z y1 y2
phi=(pi/180)*39;
delta=(pi/180)*26;
gma=18.4;
h=4;
q=20;
a=0.1;
b=0.4;
h1=1.021121318;
h2=0.744439152;
L=h+h1+h2;

beta=5;
alfa=1;
% z1=z-h;
% z2=z-h-h1;
R1=-1; 
R3=-1,%-(alfa*((z-h-h1)/h2))^0.5;
% R2=3*(beta*(1-((z-h)/(h1))))^0.5;
% R4=3*(alfa*((z-h-h1)/h2))^0.5;
kh=0.0;
kv=0;

psi=atan(kh/(1-kv));
da1=delta; pha1=phi;
l=pha1+da1;
m=pha1-psi;
n=psi+da1;

   B=(2*q*b*(b/(2*a+b)))/(gma*(h+h1)*(h+h1));
   alphacm=atan((sin(l)*sin(m)+(((sin(l)^2)*(sin(m)^2))-B*cos(n)*sin(m)*cos(l)+cos(m)*sin(l)*sin(m)*cos(l))^0.5)/((cos(m)*sin(l))-B*cos(n)))
   kg1=((tan(alphacm-phi)+tan(psi))*cos(alphacm-phi))/(tan(alphacm)*(cos(alphacm-phi-delta)))
   r1=(b/(h+h1))*(b/(2*a+b))*(tan(alphacm))
   kq1=r1*kg1
%    pgx=gma*(1-kv)*kg1*x
   pqx=(1-kv)*q*kq1
%
k4R0=(2*cos(phi-psi)^2)/(cos(phi-psi)^2+cos(psi)*cos((delta/2)+psi)...
     *(1+sqrt((sin(phi+(delta))*sin(phi-psi))/cos((delta/2)+psi)))^2)
%
delm3=-0.5*(1-R3)*delta;
k3=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R3)+cos(psi)*cos(delm3+psi)*(1-R3)*(1+sqrt((sin(phi+delm3)*sin(phi-psi))/cos(delm3+psi)))^2)
%
delm23=0.5*(3-1)*delta;
        k23=1+0.5*(3-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm23+psi)*...
            (-sqrt((sin(phi+delm23)*sin(phi-psi))/(cos(delm23+psi)))+1)^2)))-1)
       
        R2=3*(beta*(1-y1))^0.5;
        delm213=0.5*(R2-1)*delta;
        k213=1+0.5*(R2-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm213+psi)*...
            (-sqrt((sin(phi+delm213)*sin(phi-psi))/(cos(delm213+psi)))+1)^2)))-1)
    
        delm201=0.5*(1-R2)*delta;
        k201=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R2)+cos(psi)*cos(delm201+psi)...
            *(1-R2)*(1+sqrt((sin(phi+delm201)*sin(phi-psi))/cos(delm201+psi)))^2)
        
         delm43=0.5*(3-1)*delta;
         k43=1+0.5*(3-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm43+psi)*(-sqrt((sin(phi+delm43)...
             *sin(phi-psi))/(cos(delm43+psi)))+1)^2)))-1)
        
         R4=3*(alfa*y2)^0.5;
         delm413=0.5*(R4-1)*delta;
         k413=1+0.5*(R4-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm413+psi)*(-sqrt((sin(phi+delm413)...
             *sin(phi-psi))/(cos(delm413+psi)))+1)^2)))-1)
         
         delm401=0.5*(1-R4).*delta;
         k401=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R4)+cos(psi)*cos(delm401+psi)...
             *(1-R4)*(1+sqrt((sin(phi+delm401)*sin(phi-psi))/cos(delm401+psi)))^2);
         
i=0;
for z=0:0.1:L
    i=i+1;
 if(z<=h)
    p1(i)=kg1*gma*z*cos(delta)+kq1*q;
    p2(i)=0;
 elseif z<=(h+h1)
    p1(i)=kg1*gma*z*cos(delta)+kq1*q;
     for y1=0:0.01:1
      if (y1>=0 && y1<=(1-(1/beta)))
         p2(i)=-gma*(z-h)*k23*cos(delm23);
      elseif (y1>=(1-(1/beta)) && y1<=(1-(1/(9*beta))))
         p2(i)=-gma*(z-h)*eval(subs(k213*cos(delm213)));
      else (y1>=(1-(1/(9*beta))) && y1<=1)
         p2(i)=-gma*(z-h)*eval(subs(k201*cos(delm201))); 
      end
     end
 else 
    p2(i)=-k3*cos(delm3)*gma*(z-h-h1)-gma*k3*h1*cos(delta);
       for y2=0:0.01:1
      if(y2>=0 && y2<=(1/(9*alfa)))
         p1(i)=gma*(z-h-h1)*eval(subs(k401*cos(delm401)))+(gma*(h1+h)+q)*k4R0*cos(delta/2);
      elseif(y2>=(1/(9*alfa)) && y2<=(1/alfa))
         p1(i)=gma*(z-h-h1)*eval(subs(k413*cos(delm413)))+(gma*(h1+h)+q)*k4R0*cos(delta/2);
      else(y2>=(1/alfa) && y2<=1)
         p1(i)=gma*(z-h-h1)*k43*cos(delm43)+(gma*(h1+h)+q)*k4R0*cos(delta/2);
      end
       end
    
 end   
end

z=0:0.1:L;

% subplot(2,1,1);
plot(p1,z,p2,z)
grid on
set(gca, 'YDir','reverse')

Answer 1

[Raj]

"The code is taking only last values, i.e., y1 and y2=1" - No it is not taking only the last values. The for loop for both y1 and y2 runs for all values from 0 to 1 with increment of 0.1

Lets see one loop:

for y1=0:0.01:1
      if (y1>=0 && y1<=(1-(1/beta)))
         p2(i)=-gma*(z-h)*k23*cos(delm23);
      elseif (y1>=(1-(1/beta)) && y1<=(1-(1/(9*beta))))
         p2(i)=-gma*(z-h)*eval(subs(k213*cos(delm213)));
      else (y1>=(1-(1/(9*beta))) && y1<=1)
         p2(i)=-gma*(z-h)*eval(subs(k201*cos(delm201))); 
      end
end

The problem here is that you are overwriting the P2 (i) value for each value of y1 and what gets finally saved is the last computed value with y1=1. If you want to save all the P2 values for every value of y1 you will have to index it properly by using a 2D matrix.

[Akshay Pratap Singh]

Thanks @ Raj for suggestion to use 2D matrix. May you please help me in writing the same.

Your help will be highly appreciated.

Thanks in advance.

[Raj]

You can declare another indexing variable just before your 'for' loop. At the end of 'for' loop, p2 will be a 2D matrix where in a particular row, each column will be values computed for each corresponding values of y1. Something like this:

j=0;
     for y1=0:0.01:1
         j=j+1;
      if (y1>=0 && y1<=(1-(1/beta)))
         p2(i,j)=-gma*(z-h)*k23*cos(delm23);
      elseif (y1>=(1-(1/beta)) && y1<=(1-(1/(9*beta))))
         p2(i,j)=-gma*(z-h)*eval(subs(k213*cos(delm213)));
      else (y1>=(1-(1/(9*beta))) && y1<=1) % I think this should be elseif or just else without any condition. Check.
         p2(i,j)=-gma*(z-h)*eval(subs(k201*cos(delm201))); 
      end
     end
 

[Akshay Pratap Singh]

I modified the code as per your suggestion. But What I got is this:

Error using plot

Vectors must be the same length.

Please help.

[Raj]

If you have updated the code as above, p1 and p2 are 2D matrices where as Z is a 1D array. You cannot use

plot(p1,z,p2,z)

directly. Obviously the dimensions will not match. You have to modify the code to extract the data from p1 and p2 which you want to plot and then use plot command.

[Akshay Pratap Singh]

Thanks @Raj. however may you please help me by writing the code.

Thanks in advance.

[Raj]

Depends on which value of y1/y2 you want to plot p2/p1 versus z. Its quite straight forward actually. p1/p2 are 2D matrices where rows are values with varying 'z' and columns are values with varying y2/y1 respectively. Just extract the relevant data and plot against 'z'.