How do I use fsolve n*m times in a 3-dimensional array without messing up my dimensions?

I am searching for solutions to a system of three equations, three unknowns. The system is called eqs and is a function of a(1), a(2), a(3), and the parameter S. For n values of the parameter S, I want to solve eqs m times with m random guesses. Here is the function I am trying to solve m*n times

    function F=eqs(a,S) 
F(1)=(S-2*(a(3))^2)*(pi-acos(a(1)/a(3))-acos(a(2)/a(3))-acos(((1-2*a(1))^2+...
    (1-2*a(2))^2)^0.5/((2*a(3)))))-a(1)*((a(3))^2-(a(1))^2)^0.5-a(2)*...
    ((a(3))^2-(a(2))^2)^0.5-1/4*((1-2*a(1))^2+(1-2*a(2))^2)^0.5*...
    (4*(a(3))^2-(1-2*a(1))^2-(1-2*a(2))^2 )^0.5;
F(2)=(4*a(2)*(a(3))^2+a(2)-4*(a(3))^2)/(4*a(2)-4*(a(3))^2-1)-a(1);
F(3)=(1-2*a(2))^2/(1-2*a(1))^2*(4*(a(3))^2-4*(a(1))^2 )-4*(a(1))^2+...
    (1-2*a(1))^2+(1-2*a(2))^2;
end

Here is the code that is not working. I have some problem with my matrix dimensions on the solutions line, but I cannot figure it out. please help.

    %Preallocations
n=5;
m=10;
draw=100;
S=zeros(n,1);
guess=zeros(n,m,3);
b=zeros(n,m,1);
c=zeros(n,m,1);
d=zeros(n,m,1);
solutions=zeros(n,m,3);
for t=1:n;
    S(t)=(t+n/(3+2^1.5))/n; %this is the parameter that varies. n values.
      for k=1:m;
    %for every value of S I want to fsolve m times with m random guesses
    %generating the guess matrix
    b(t,k)=randsample(draw,1);
    c(t,k)=randsample(draw,1);
    d(t,k)=randsample(draw,1); 
    guess(t,k,1)=b(t)/(2*draw+1);
    guess(t,k,2)=c(t)/(2*draw+1);
    guess(t,k,3)=d(t)*1.2/(2*draw+1);
    %solving the off-diagonal system
    system = @(a)eqs(a,S);
    solutions(t,k,:)=fsolve(system,guess(t,k,:));%=[a(1) a(2) a(3)]
  end
end

Answer 1

[Star Strider]

You need to subscript ‘S’ in your function call to pass it to your function as a scalar:

system = @(a)eqs(a,S(t));

and your code works. Otherwise, it is a vector and it will throw an error. I leave it to you to determine if it gives the results you want.

[Matt J]

I recommend using fminsearch instead,

solutions(t,k,:) = fminsearch( @(a) norm( eqs(a,S(t)) ),  guess(t,k,:));

since eqs() is not differentiable, and it is not obvious to me how to rewrite it so that it is.

[Matt]

I'll have to look into that. It is differentiable, but it is UGLY. And thank you Star Strider! That did it!

[Star Strider]

As always, my pleasure!

[Matt J]

It is differentiable, but it is UGLY.

I see expressions like this,

((a(3))^2-(a(1))^2)^0.5

which are non-differentiable wherever a(3)=a(1). You would have to be sure that fsolve will not look for a solution in that region.