Without using drawnow command i want plot from below code

function run_LE_FO_p1(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p_min,p_max,n);
figure();
hold on;
ne=3;
ext_fcn=@LE_RF_p1;
t_start=0;h_norm=0.02;t_end=10;
x_start=[0.1 0.1 0.1];h=0.02;q=0.998;
p_max=1.3;
p_min=1.1;
n=800;
p_step=(p_max-p_min)/n
p1=p_min;
while p1<=p_max
    lp=FO_Lyapunov_p1(ne,ext_fcn,0,0.02,10,[0.1 0.1 0.1],0.002,0.998,p1)
    p1=p1+p_step;
    plot(p1,lp,'.')
    drawnow();
end
function LE=FO_Lyapunov_p1(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p1);
% Memory allocation
x=zeros(ne*(ne+1),1); 
x0=x;
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);
% Initial values
x(1:ne)=x_start;
i=1;
while i<=ne
    x((ne+1)*i)=1.0;
    i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
    % Solutuion of extended ODE system of FO using FDE12 routine
    [T,Y] = FDE12(q,ext_fcn,t,t+h_norm,x,h,p1); 
    t=t+h_norm;
    Y=transpose(Y);
    x=Y(size(Y,1),:); %solution at t+h_norm
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x0(ne*i+j)=x(ne*j+i);
            j=j+1;
        end;
        i=i+1;
    end;
    %       construct new orthonormal basis by gram-schmidt
    zn(1)=0.0;
    j=1;
    while j<=ne
        zn(1)=zn(1)+x0(ne*j+1)^2; 
        j=j+1;
    end;
    zn(1)=sqrt(zn(1));
    j=1;
    while j<=ne
        x0(ne*j+1)=x0(ne*j+1)/zn(1); 
        j=j+1;
    end
    j=2;
    while j<=ne
        k=1;
        while k<=j-1
            gsc(k)=0.0;
            l=1;
            while l<=ne;
                gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
                l=l+1;
            end
            k=k+1;
        end
        k=1;
        while k<=ne
            l=1;
            while l<=j-1
                x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);
                l=l+1;
            end
            k=k+1;
        end;
        zn(j)=0.0;
        k=1;
        while k<=ne
            zn(j)=zn(j)+x0(ne*k+j)^2;
            k=k+1;
        end
        zn(j)=sqrt(zn(j));
        k=1;
        while k<=ne
            x0(ne*k+j)=x0(ne*k+j)/zn(j);
            k=k+1;
        end
        j=j+1;
    end
    %       update running vector magnitudes
    k=1;
    while k<=ne;
        c(k)=c(k)+log(zn(k));
        k=k+1;
    end;
    %       normalize exponent        
    k=1;
    while k<=ne
        LE(k)=c(k)/(t-t_start);
        k=k+1;
    end
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x(ne*j+i)=x0(ne*i+j);
            j=j+1;
        end
        i=i+1;
    end;
    x=transpose(x);
    it=it+1;
end
function f=LE_RF_p1(t,x,p1)
%p is the parameter
f=zeros(9,1);
X= [x(4) x(7) x(10);
    x(5) x(8) x(11);
    x(6) x(9) x(12)];
%RF equations
f(1)=x(2).*(x(3)-1+x(1).*x(1))+0.1*x(1);
f(2)=x(1).*(3*x(3)+1-x(1).*x(1))+0.1*x(2);
f(3)=-2.*x(3).*(p1+x(1).*x(2));
%Jacobian matrix
J=[2*x(1).*x(2)+0.1, x(1).*x(1)+x(3)-1, x(2);
    -3*x(1).*x(1)+3*x(3)+1,0.1,3*x(1);
    -2*x(2).*x(3),-2*x(1).*x(3),-2*(x(1).*x(2)+p1)];
f(4:12)=J*X; % To be modified if ne>3

I am getting output as dotted and non fluctuated plot

i want plot like below attachment

please give me the suggestion