Нужно дополнить код решением, используя классический метод Ньютона

iBis=0;%constants
imNM=0;  
e=0.0000001;
lb=-10;
rb=-4;
pkx1=-10;
pk=-4;
 
p=-10:0.1:-4;
%fzero
f=@(p) (0.001*p.^5)+(0.02*p.^4)+(0.5*p.^3)+(4*p.^2)+10*p+40; 
y=f(p);
figure;
plot(p,y);
grid;
p1=[-10 -4];
opts=optimset('Display','iter','Tolx',1e-5);%TolX – final p tolerance
disp('Information about function fzero'); 
[p, fval]=fzero(f,p1,opts);
% p01=fzero(f, p1);
hold on;
plot(p,fval,'ko');
 fprintf('Root (fzero): %.8f\n', p); %formatted output
%Bisection Method
while (abs(rb-lb)>e) 
 iBis=iBis+1; %counter
 mid=(rb+lb)/2;
 if(f(lb)*f(mid)<0)
     rb=mid;
 elseif (f(mid)*f(rb)<0)
     lb=mid;
 else
     break;
 end
end
fprintf('Root (BisectionMethod): %.8f\n', mid);
fprintf('Number of iterations (BisectionMethod): %d\n', iBis);
%Modified Newton's Method (mNM)
%finding differential 
 syms x;           %create symbolic variable
 y=(0.001*x.^5)+(0.02*x.^4)+(0.5*x.^3)+(4*x.^2)+10*x+40; %my function, but p==x
 z1=diff(y);       %finding differential of my function
                   %it returns the equation in symbolic form
 z2=diff(z1);
 z2=subs(z2,x,pkx1);
 z2=double(z2);
 z1=subs(z1,x,pkx1);%replace x with pkx1(right border) 
 z1=double(z1);     %convert to numerical form
if (f(pkx1)*z2>0)
 while(abs(pkx1-pk)>e)
     imNM=imNM+1; %counter
     pk=pkx1;     %write the previous value
     pkx1=pk-f(pk)/z1; %count next value
     
     if (f(pkx1)==0) %if funсtion = 0 then 
        break;       %we stop the cycle 
     end
 
 end
 fprintf('Root (mNM): %.8f\n',pkx1);%and display the data
 fprintf('Number of iterations (mNM): %d\n',imNM);
 else
     return;
end