MATLAB: FSOLVE requires all values returned by functions to be of data type double.

fsolveMATLAB

while running this code i am getting "FSOLVE requires all values returned by functions to be of data type double" error.
Description:
for a different value of "a" i want to calculate b & c for plotting a graph between 'sigmaxx' vs 'str.
str= a (from 0.5 to 1.5)
I understand that i have to to change the symbolic function to numeric but i am not able to do using double.
what are the possible & efficient way to change the variable from symbolic to numeric in the following code.
Thanks in advance.
Code:
clc
clear all
close all
syms b c
a=linspace(0.5,1.5,100);
nu=0.5;
gamm=2e6;
k=0;
I=eye(3,3);
for p=1:length(a)
    F=[a(p) k 0;0 b 0;0 0 c];
    B=F*F'; 
    Binv=inv(B);
    J1=trace(B);
    J2=trace(Binv);
    J3=det(F);
    sigma(:,:,p)=-nu*J2*J3*I+(((nu*(J1-1)/J3) + (J1-3/(2*J3))*gamm))*B +nu*J3*Binv;
    sigmaxx(p)=sigma(1,1,p);
    sigmayy(p)=sigma(2,2,p);
    sigmazz(p)=sigma(3,3,p);
    g=@(b,c) sigmayy(p);
    h=@(b,c) sigmazz(p);
    f=@(w) [g(w(1),w(2));h(w(1),w(2))];
    w0=[0.6683,0.6683];                   % value of b & c when a=1.3
    options=optimset('MaxIter',1e5,'TolFun',1e-10);
    wroot=fsolve(f,w0,options);
    b=wroot(:,1)
    c=wroot(:,2)
end
str=a;
plot(str,sigmaxx)

Best Answer

I don't see any need for syms.
a=linspace(0.5,1.5,100);
nu=0.5;
gamm=2e6;
k=0;
I=eye(3,3);
w0=[0.6683,0.6683];                   % value of b & c when a=1.3
options=optimset('MaxIter',1e5,'TolFun',1e-10,'Display','none');
[b,p]=deal(nan(1,length(a)));
for p=1:length(a)
    f=@(w) func(w,a(p),gamm,nu,k,I);
    wroot=fsolve(f,w0,options);
    
    b(p)=wroot(1);
    c(p)=wroot(2);
    
end
whos b c
  Name      Size             Bytes  Class     Attributes

  b         1x100              800  double              
  c         1x100              800  double              
function out=func(w,a,gamm,nu,k,I)
    b=w(1); c=w(2);
    F=[a k 0;0 b 0;0 0 c];
    B=F*F'; 
    Binv=inv(B);
    J1=trace(B);
    J2=trace(Binv);
    J3=det(F);
    sigma=-nu*J2*J3*I+(((nu*(J1-1)/J3) + (J1-3/(2*J3))*gamm))*B +nu*J3*Binv;
    out(1)=sigma(2,2);
    out(2)=sigma(3,3);
end
Related Question