I want to plot phase diagram .But I reach here maintion the code below please help ,me the remaining part

clc;clear;close all;

r=0.1;

k=50;

a=0.01;

e=0.5;

m=0.05;

s=0.1;

w=0.1;

b=0.001;

q=0.03;

F=1.;

M=50;

X=0:1:50;J=[0 50 0 20];

Pi_0=(w/(w+b*M))*(F/s);Pi_1=((w+b*M)/w)*(F/s);

D1=Pi_0+0*X;D2=Pi_1+0*X;

plot(D1,X,'k--')

axis(J)

hold on

plot(D2,X,'k--')

X1=(Pi_0/99);X2=(Pi_1-Pi_0)/99;X3=(50-Pi_1)/99;Y2=20/99;

P1=0:X1:Pi_0;H1=0:Y2:20;P2=Pi_0:X2:Pi_1;H2=0:Y2:20;P3=Pi_1:X3:50;H3=0:Y2:20;

u_1=0;u_2=((s.*P2)/(s.*P2+F)+(w*(s.*P2-F))/(b*M*(s.*P2+F)));u_3=1;

f = @(P1,H1) (r.*P1).*(1-P1/k)-(a.*P1.*H1)./(1+q*u_1*M);

fimplicit(f,[0 Pi_0 0 20],'g')

hold on

f1 = @(P1,H1) (e*a.*P1.*H1)./(1+q*u_1*M)- m.*H1;

fimplicit(f1,[0 Pi_0 0 20],'g')

hold on

%%

[P1, H1] = meshgrid(0, 0:0.01:Pi_0);

% dP1dt = dP1(P1, H1);

% dH1dt = dH1(P1, H1);

quiver(P1, H1, dP1dt, dH1dt, 'color', 'b');

f = @(P2,H2) (r.*P2).*(1-P2/k)-(a.*P2.*H2)./(1+q*((s.*P2)/(s.*P2+F)+(w*(s.*P2-F))/(b*M*(s.*P2+F)))*M);

fimplicit(f,[Pi_0 Pi_1 0 20],'r')

hold on

f1 = @(P2,H2) (e*a.*P2.*H2)./(1+q*((s.*P2)/(s.*P2+F)+(w*(s.*P2-F))/(b*M*(s.*P2+F)))*M)- m.*H2;

fimplicit(f1,[Pi_0 Pi_1 0 20],'r')

hold on

%%

[P2, H2] = meshgrid(0:0.01:Pi_0, 0:0.01:Pi_1);

% dP2dt = dP2(P2, H2);

% dH2dt = dH2(P2, H2);

quiver(P2, H2, dP2dt, dH2dt, 'color', 'b');