Muller Method Code in Matlab (Easy to Understand)

syms x a b  c;
f=16*x^4 - 40*x^3 + 5*x^2 + 20*x + 6;
syms p0 ;syms  p1 ; syms p2 ;
p0=0.5;p1=-0.5;p2=0;
iterate=10;
error=[];
p=[];

for i=1:iterate
c=subs(f,p2); %value of cd
a=(((p1-p2)*(subs(f,p0)-subs(f,p2)))-((p0-p2)*(subs(f,p1)-subs(f,p2))))/((p0-p1)*(p0-p2)*(p1-p2));
b=((((p0-p2)^2)*(subs(f,p1)-subs(f,p2)))-(((p1-p2)^2)*(subs(f,p0)-subs(f,p2))))/((p0-p1)*(p0-p2)*(p1-p2));
disp('p3');
p3=p2-(2*c/(b+(sign(b)*(sqrt(b^2-4*a*c)))));
disp(p3);

error(i)=p3-p2;
p0=p1;
p1=p2;

p2=p3;