This is Broyden.m I found. I have a system of 3 functions and 6 unknowns, I think this function is meant to solve 1 unknown per equation or am I not coding my eqns() function corrrectly?(in comments)

function [x, ithist] = broyden(f,x0,opt,bounds)

% Solve system of nonlinear equations by Broyden's method

% Broyden's method is a quasi Newton method which updates an approximate

% Jacobian at each new Newton step.

%

% Usage:

% [x, ithist] = broyden( f, x0, opt, bounds)

%

% Input parameters:

% f name of a matlab function that evaluates F. (Required)

% x0 iteration starting point (Required)

% opt Options struct. Fields: (Optional)

% tolfun Function tolerance. (Default tolfun = 1.e-10)

% Iterations end if ||F(x)||_2 < tolfun

% tolx Step length tolerance (Default tolx = 1.e-8)

% Newton's method stops if ||step||_2 < tolx,

% where step is the Newton step.

% maxiter maximum number of iterations (Default maxiter = 50)

% bounds Optional bounds on x Default: no bounds.

% bounds(:,1) <= x < bounds(:,2)

%

% Output parameters:

% x approximation of the solution.

% ithist iteration history struct. Fields: x, f(x), norm(f)

%

% Example:

% f = @(x) [x(1)^2 + x(2)^2 - 4; exp(x(1)) + x(2) - 1];

% x = broyden(f,[1;1])

% The other root has x(1) > 0, x(2) < 0. This may be found by selecting

% another inital x, or by using constraints. E.g.:

% x = broyden(f,[1;1],[],[0 2.5;-2.5 0])

%

% See also: fsolve (in optimization toolbox)

% Author: Are Mjaavatten, Telemark University College, Norway

%

% The update formula for the Jacobian J is taken from

% the broyden function written by Matthias Heinkenschloss:

% http://read.pudn.com/downloads158/sourcecode/math/705810/sysnlineq/broyden.m__.htm

%

% Version 1: 2015-12-29

% Version 1.1: 2016-11-29 Added iteration history

optfields = {'maxiter','tolfun','tolx'};

defaults = {50,1e-10,1e-8,[]};

if nargin < 3

opt = [];

end

for i = 1:3

if isfield(opt,optfields{i})

if isempty(opt.(optfields{i}))

opt.(optfields{i}) = defaults{i};

end

else

opt.(optfields{i}) = defaults{i};

end

end

x = x0(:);

it = 0;

F = feval(f, x);

if ~(size(x,1) == size(F,1))

error('f must return a column vector of the same size as x0')

end

normf = norm(F);

J = jacobi(f,F,x); % Intial Jacobian matrix

if nargout > 1

ithist.x = [x(:)';zeros(opt.maxiter,length(x))];

ithist.f = [F(:)';zeros(opt.maxiter,length(x))];

ithist.normf = [normf;zeros(opt.maxiter,1)];

end

normdx = 2*opt.tolx;

while( it < opt.maxiter+1 && normdx > opt.tolx && normf > opt.tolfun)

if rcond(J) < 1e-15

error('Singular jacobian at iteration %d\n',it)

end

dx = -J\F;

normdx = norm(dx);

if nargin > 3 % variable bounds are supplied

% make sure x stays within bounds

for j = 1:20

jl = find(x+dx<bounds(:,1));

dx(jl) = dx(jl)/2;

ju = find(x+dx>bounds(:,2));

dx(ju) = dx(ju)/2;

if isempty(jl) && isempty(ju)

break

end

end

end

x = x+dx;

it = it+1;

F = feval(f, x);

normf = norm(F);

J = J + F*dx'/(dx'*dx);

if nargout > 1

ithist.x(it+1,:) = x(:)';

ithist.f(it+1,:) = F(:)';

ithist.normf(it+1,:) = normf;

end

end

% Check if the iterations converged and issue warning if needed

if it >= opt.maxiter && norm(F) > opt.tolfun

warning('No convergence in %d iterations.\n',it+1)

elseif normf>opt.tolfun

warning('Newton step < %g, but function norm > %g\n',...

opt.tolx,opt.tolfun)

elseif normdx>opt.tolx

warning('Function norm < %g, but newton step norm > %g\n',...

opt.tolfun,opt.tolx)

end

if nargout > 1

ithist.x(it+2:end,:) = [];

ithist.f(it+2:end,:) = [];

ithist.normf(it+2:end) = [];

end

end

function J = jacobi(f,y0,x)

% Quick and dirty numerical Jacobian for function f at x

% y0: f(x);

delta = 1e-6*(max(1,sqrt(norm(x))));

n = length(y0);

m = length(x);

J = zeros(n,m);

for i = 1:m

dx = zeros(m,1);

dx(i) = delta/2;

J(:,i) = (feval(f,x+dx)-feval(f,x-dx))/delta;

end

end

[Timothy Nsubuga]

%*************************************************************************

%%This fucntion will store transducer coordinates and differences of

%%arrival times between each transducer and transducer 1,

%the arrival times are aquired via graphData() and is an input broyden() to solve for the

% signal location and travel time to transducer 1, adding the difference in

% arrival time of each transducer to this travel time yields the travel

% time

%%

%**************************************************************************

function [fun] = eqns(u)%unknown values are input parameters

%the pop rock falls vertically and so xs = zs = 0 and we assume the center

%of the borehole is at y = 0

xs = 0;

zs = 0;

%Transducer 1 and 2 locations in meters based on 3D print sample with a

%fracture 1/3 of the sample height

xt1 = -0.05;

yt1 = 0.015;

zt1 = -0.015;

xt2 = -0.05;

yt2 = 0.015;

zt2 = 0.015;

xt3 = -0.05;

yt3 = -.005;

zt3 = -0.015;

v = 2500;% wave speed in material determined by AST and distance formula

%differnce in arrival time between tranducer 1 and transducer 2;

% acquire locations of first peaks in transducers 1 and 2

File_names = dir('*.csv'); %attain every waveform in the event folder

%get files of interest, in this case those of transducer 1 and 2

file1 = File_names(3).name;

file2 = File_names(4).name;

file2 = File_names(5).name;

plus1 = 1; %offset y values

plus2 = 2;

plus3 = 3;

hold on

% %prof requested the dts so I need all the first peak locations

% numFiles = length(File_names);

% incr = 1:numFiles;

% for i = 1:numFiles %call graphData() for each waveform

% file = File_names(i).name;

% hold on

% plus = incr(i);%needed to offset graphs for plotting

% [p,l] = graphData(file,plus);

%

%

% end

[p1,l1] = graphData(file1,plus1);

% (ORIGINAL CODE IS COMMENTED OUT FOR NOW TO GET ALL

% DTS FROM EVENT 1)

% L1 = l1;

[p2,l2] = graphData(file2,plus2);

[p3,l3] = graphData(file2,plus2);

% L2 = l2;

% graphData() gives me locations on the first peak but I only require one of

% them

% maxL1 = max(L1);

%

% maxL2 = max(L2);

%multiply each location by the sample interval, 1*10^-7 to get the arrival

%times for transducers 1 and 2

art1 = l1*1e-7;

art2 = l2*1e-7;

art3 = l3*1e-7;

%

% tran10= 1038*1e-7;

%

% tran13 = 1019*1e-7;

%

% tran1 = 1047*1e-7;

%

% tran2 = 1054*1e-7;

%

% tran3 = 1038*1e-7;

%

% tran5 = 1048*1e-7;

%

% tran6 = 1042*1e-7;

%

% tran8= 1022*1e-7;

%

% tran9 = 1034*1e-7;

%

% dt10 = tran10 - tran1;

% dt13 = tran13 -tran1;

% dt1 = tran1- tran1;

% dt2 = tran2- tran1;

% dt3 = tran3- tran1;

% dt5 = tran5- tran1;

% dt6 = tran6- tran1;

% dt8 = tran8- tran1;

% dt9 = tran9 - tran1;

%deltat1 is clearly zero

deltat2 = art2 - art1;

deltat3 = art3 - art1;

% deltat2 = 7e-7;

ys = u(1); %source location and travel time from source to transducer 1

% are unkown and will be guessed initally with values of

% [0.5;2e-5] as input argumnents to eqns()

ts = u(2);

%non-linear system of equations, the first is for transducer 1 and the

%second is for transducer 2

fun(1,1) = (sqrt((xt1- xs)^2 + (yt1 - ys)^2 + (zt1 - zs)^2)/v )- (ts);

fun(2,1) = (sqrt((xt2- xs)^2 + (yt2 - ys)^2 + (zt2 - zs)^2)/v) - (ts - deltat2);

fun(3,1) = (sqrt((xt3- xs)^2 + (yt3 - ys)^2 + (zt3 - zs)^2)/v) - (ts - deltat3);

end

%end of eqns()*************************************************************

[Walter Roberson]

You are right that the code cannot have a different number of variables than it has return values from the function.

It might be interesting to return 3 extra 0 from the function.

[Timothy Nsubuga]

guess = [0.5;2e-5];

>> f = @eqns;

>> x = broyden(f,guess);

Here is my call, each equation has ts and ys, and I want to solve for these 2 variables. guess is the corresoponding values I beleive ys and ts are, any tips to get this to work?