function [xList, functionValueList] = Mullers(f, p0, p1, p2, errorAllowed)
%
%Mullers(f, p0, p1, p2, errorAllowed)
%This function finds the root of a function using Muller's Method
%

	%Make sure the number of arguments is correct
	if(nargin ~= 5)
		error('That is an incorrect number of arguments')
	end
	%A few necesary things before we get started
	pkg load symbolic;
	warning('off','OctSymPy:sym:rationalapprox');
	%Constant variables
	maxIt = 50;
	%Variables
	h1 = p1 - p0;
	h2 = p2 - p1;
	s1 = (double(subs(f,p1)) - double(subs(f,p0))) / h1;
	s2 = (double(subs(f,p2)) - double(subs(f,p2))) / h2;
	d = (s2 - s1) / (h2 + h1);
	cnt = 2;

	%Loop until you reach the maximum number of itterations
	%If you find an answer within error it will return
	while(cnt < maxIt)
		b = s2 + h2 * d;
		D = (b^2 - (4 * double(subs(f,p2)) * d))^(1/2);
		if(abs(b - D) < abs(b + D))
			E = b + D;
		else
			E = b - D;
		end
		h = (02 * double(subs(f,p2)))/E;
		p = p2 + h;
		xList(end+1) = p;
		functionValueList(end+1) = double(subs(f,p));
		%Exit the function if you get an answer within error
		if(abs(h) < errorAllowed)
			return;
		else
			p0 = p1;
			p1 = p2;
			p2 = p;
			h1 = p1 - p0;
			h2 = p2 - p1;
			s1 = (double(subs(f,p1)) - double(subs(f,p0)))/h1;
			s2 = (double(subs(f,p2)) - double(subs(f,p1)))/h2;
			d = (s2 - s1)/(h2 + h1);
			++cnt;
		end
	end
end