Plot dynamics of an ODE

I have an ODE and I need to plot with respect to and varying p. Following is the sample plot required.

I managed to draw one plot for and following is the sample code.

[xi,x] = meshgrid(0:0.001:1,0:0.002:1);
p =1
f2 = @(x,a,b) a*b*(x)./(b-2+x));
f1 = @(x,a,b) a*b*(1-x)./(2*b-x);

rs = @(x1,x2,a,b,p) x1.*f1(x,a,b)  -   p.*(1-x1).*f2(x,a,b);

figure
contour(x1,2,rs(x1,x2,1,2,1),[0 0]);

Next I added a loop for p as follows:

np = 10;
phi = linspace(0,1,np);

for p = 1:np
    f2 = @(x,a,b) a*b*(x)./(b-2+x));
    f1 = @(x,a,b) a*b*(1-x)./(2*b-x);

    rs = @(x1,x2,a,b,p) x1.*f1(x,a,b)  -  p.*(1-x1).*f2(x,a,b);
    R = rs(x1,x2,1,2, p(p));
    contour(x1,x2,R);
    hold on
end
hold off 

But, instead of the required ourput of 10 plots in one graph, I get many plots (>10) in my output. Can someone please help me with this?

[Cris LaPierre]

You have not copied your code correctly here, as both examples do not run.

Answer 1

[Cris LaPierre]

I'm assuming your goal is to plot a line for when the function rs is equal to zero for various values of p. You could do that using contour, but I don't think that is the best way. However, you could still use contour to get the (x,y) coordinates of the 0 contour line for each, and then plot those.

np = 10;
[x1,x2] = meshgrid(0:0.001:1,0:0.002:1);

for p = 1:np
    f2 = @(x,a,b) a*b*(x)./(b-2+x);
    f1 = @(x,a,b) a*b*(1-x)./(2*b-x);

    rs = @(xi,x,a,b,p) xi.*f1(x,a,b)  -  p.*(1-xi).*f2(x,a,b);
    R = rs(x1,x2,1,2, p);
    figure(1)
    M=contour(x1,x2,R,[0 0]);

    figure(2)
    plot(M(1,2:end),M(2,2:end),'DisplayName',"p="+p)
    axis([0 1 0 1])
    hold on
end

hold off 
legend(Location="northwest")

[UserCJ]

@Cris LaPierre Thank you very much for your answer. Can you also tell me the best way to plot a line for when the function rs is equal to zero for various values of p? And I'm thinking of looping over x1 and x2 instead of using a meshgrid as shown below;

N1 = 10;
N2 = 10;
x1 = linspace(0,1,N1);%x1
x2 = linspace(0,1,N2);%x2

nphi = 20;
phi = linspace(0,1,nphi);
a=1;
b=3; 
dxdt_A = zeros(N1,N2);

f1 = @(x,a,b) a*b*(x)/(b-2+x);
f2 = @(x,a,b) a*b*(1-x)/(2*b-x);

% 
for p = 1:nphi
    for i = 1:N1
        for j = 1:N2
            dxdt_A(i,j) = f1(x2(j),x1(i),b,a) - phi(p) * f2(x2(j),x1(i),b,a);

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

        end
    end
  M=contour(x1,x2,dxdt_A,[0,0]);
  figure(1)
  plot(M(1,2:end),M(2,2:end),'DisplayName',"phi="+p,'Linewidth',2)
  hold on 
end
hold off

But this doesn't give me a good plot! Can you please help me with this?

[UserCJ]

@Cris LaPierre It was an error in parantheses. I was experimenting on transforming variables to a different domain, so forgot to remove one bracket. sorry about that! However, I still do not get the same plot.

[Cris LaPierre]

I would keep the meshgrid approach. That will be faster than looping.

I'm not sure of a best way, but one way to do this without contour is to look for the min(abs(R)) in each row. This won't be the exact zero point, but will be close enough for the visualization.

np = 10;
X = 0.001:0.001:1;
Y = 0.002:0.002:1;
[x1,x2] = meshgrid(X,Y);

for p = 1:np
    f2 = @(x,a,b) a*b*(x)./(b-2+x);
    f1 = @(x,a,b) a*b*(1-x)./(2*b-x);

    rs = @(xi,x,a,b,p) xi.*f1(x,a,b)  -  p.*(1-xi).*f2(x,a,b);
    R = rs(x1,x2,1,2, p);
    
    [xval,xind] = min(abs(R),[],2);
    X_z0(:,p) = X(xind);
end

plot(X_z0,Y)
axis([0 1 0 1])
legend(Location="northwest")

[Cris LaPierre]

You may find this answer helpful as well

https://www.mathworks.com/matlabcentral/answers/480655-how-to-find-the-x-y-combinations-that-make-z-0