Help with my secant method.

Here's my code.

function [it_count,root,xn] = secant(fcn,x0,x1,error_bd,max_iterate)
%
% function secant(fcn,x0,x1,error_bd,max_iterate)
%
% This implements the secant method for solving an 
% equation f(x) = 0. The function f=fcn is a string.
%
% The parameter error_bd is used in the error test for the 
% accuracy of each iterate.  The parameter max_iterate is an 
% upper limit on the number of iterates to be computed. Two 
% initial guesses, x0 and x1, must also be given.
%
% For the given function f(x), an example of a calling sequence 
% might be the following:
%    root = secant('x^2-1',x0,x1,1.0E-12,10)
%
% The program optionally prints the iteration values
%      iterate_number, x, f(x), error
% The value of x is the most current initial guess, called
% previous_iterate here, and it is updated with each iteration. 
% The value of error is 
%    error = newly_computed_iterate - previous_iterate
% and it is an estimated error for previous_iterate.

error = 1;
f=inline('x^5 - x^4 + x -1');
fx0 = feval(f,0.5); %feval --> fzero
disp("fx0") %**********
disp(fx0)   %**********
it_count = 0;
%format short e
%iteration = [it_count x0 fx0]
xsave=zeros(100,1); xn=0; % for saving each iteration
xsave(1)=0.5;
xsave(2)=2;

while abs(error) > 1.0E-6 & it_count <= 100
    fx1 = feval(f,2);
    if fx1 - fx0 == 0
        disp('f(x1) = f(x0); Division by zero;  Stop')
	root=NaN;
        return
    end
    it_count = it_count + 1;
    x2 = 3 - fx1*(2-0.5)/(fx1-fx0);
    xsave(it_count+2)=x2;
    error = x2 - 2;
    %Internal print of secant method:
    %format short e
    iteration = [it_count 2 fx1 error]
    x0 = 2;
    x1 = x2;
    fx0 = fx1;
%     disp("fx0_1") %**********
%     disp(fx0)   %**********
%     disp("fx1_1")
%     disp(fx1)   %**********
%     disp("x1_1")
%     disp(x1)   %**********
%     disp("x2_1")
%     disp(x2)    %**********
end

root = x2;
xn=xsave(1:it_count+2);
disp("fx0_1") %**********
disp(fx0)   %**********

if it_count > 100
  disp('The number of iterates calculated exceeded') 
  disp('max_iterate.  An accurate root was not')
  disp('calculated.')
end

% % disp("test_last")
% % disp(x0)
% % disp(x1)
% % disp(iteration)
% % disp(fx0)

%format long
%root
%format short e
%error
%format short
%it_count

Output:

secant

fx0

-0.5312

iteration =

1.0000 2.0000 17.0000 -0.4545

f(x1) = f(x0); Division by zero; Stop

ans =

1

The result shws only 1 iteration.

How do I fix this?

Answer 1

[darova]

Some suggestions:

x0 = 0.5;
x1 = 2;
fx0 = feval(f,x0);
while abs(error) > 1.0E-6 & it_count <= 100
%     fx1 = feval(f,2);                     % why 2 is a constant? should be x1
    if fx1 - fx0 == 0
        disp('f(x1) = f(x0); Division by zero;  Stop')
    root=NaN;
        return
    end
    it_count = it_count + 1;
    x2 = 3 - fx1*(2-0.5)/(fx1-fx0);        % why 3, 2 and 0,5 are constants?
                                           % i'd try: x2 = x1 - fx1*(x1-x0)/(fx1-fx0)
    xsave(it_count+2)=x2;
    error = x2 - 2;                        % don't know about this, maybe error = fx2 ?
    %Internal print of secant method:
    %format short e
    iteration = [it_count 2 fx1 error]
    
%   I'd try:
%         fx2 = feval(f,x2);
%         if fx2 < 0
%             x0 = x2;
%             fx0 = fx2;
%         else
%             x1 = x2;
%             fx1 = fx2;
%         end

%     x0 = 2;                                % why always 2?
%     x1 = x2;
%     fx0 = fx1;
end

Also

f=inline('x^5 - x^4 + x -1');
fx0 = feval(f,0.5); %feval --> fzero

f = @(x) x.^5 - x.^4 + x -1;            % it's better and faster
fxx = f(x0);

Answer 2

[Nattandige Fernando]

[Walter Roberson]

No, that does not appear to be an Answer to the Question that was asked.