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

Hello. this is my code.it solves the first function and gives the values of x1(i,1).

but it does not go to solving the function of mysubfun1. and stops.

what should i do?

function x1 = supercriticalahmad

for i=1:2

x1(i,1) = fsolve(@mysubfun,0.1);

p2x1(i,1)=sqrt(1-(4.*x1(i,1).^2));

p1x1(i,1)=2.*asin(2.*x1(i,1))+4.*x1(i,1).*sqrt(1-(4.*x1(i,1).^2));

F1(i,1)=(p2x1(i,1)/p2xc(i,1))^(1/2).*(p1xc(i,1)/p1x1(i,1))^(3/2);

end

function xe=hoda

xe(i,1) = fsolve(@mysubfun1,0.1)

EDR(i,1)=xe(i,1)./xc(i,1)

function ff= mysubfun1(xe)

p1xe(i,1)=2.*asin(2.*xe)+4.*xe.*sqrt(1-(4.*xe.^2));

b(i,1)=1/(2.*xc(i,1).*(1+(p1xc(i,1)/(16.*xc(i,1).*p2xc(i,1))))); %b(i,1) haman R/H mib(i,1)ashad.

a0(i,1)=1;

a1(i,1)=(1/2).*b(i,1);

a2(i,1)=(1/8).*((b(i,1).^2)+4);

a3(i,1)=(1/16).*((b(i,1).^3)-(4.*b(i,1)));

a4(i,1)=(5/128).*((b(i,1).^4)-(8/5).*b(i,1).^2+(16/5));

a5(i,1)=(7/256).*(b(i,1).^5-(8/7).*b(i,1).^3-(16/7).*b(i,1));

a6(i,1)=(21/1024).*((b(i,1).^6)-(20/21).*b(i,1).^4-(16/21).*b(i,1).^2+(64/21));

a7(i,1)=(231/14336).*((b(i,1).^7)-(28/33).*b(i,1).^5-(112/231).*b(i,1).^3-(448/231).*b(i,1));

p3x1(i,1)=((a0(i,1).*xc(i,1))-(a1(i,1).*(xc(i,1).^2))-((4/3).*a2(i,1).*(xc(i,1).^3))...

-((8/4).*a3(i,1).*(xc(i,1).^4))-((16/5).*a4(i,1).*(xc(i,1).^5))-((32/6).*a5(i,1).*(xc(i,1).^6))...

-((64/7).*a6(i,1).*(xc(i,1).^7))-((128/8).*a7(i,1).*(xc(i,1).^8)));

ff=(32.*sqrt(2).*p1xe(i,1).*p3x1(i,1).*((x1(i,1)+(F1(i,1).^2.*p1x1(i,1))...

/(16.*p2x1(i,1))).^(1/2)).*(p2x1(i,1).^(1/2)./p1x1(i,1).^(5/2)))-F1(i,1);

end

end

disp('xc=')

disp(xc)

disp('Qstar=')

disp(Qstar)

disp('sc=')

disp(sc)

disp('x1=')

disp(x1)

disp('F1=')

disp(F1)

%disp('xe=')

%disp(xe)

%disp('EDR=')

%disp(EDR)

function f = mysubfun(x1)

Q=[0.000788

0.000866];

s=0.022;

n=0.012;

D=0.128;

xc=[0.123

0.131];

p1xc(i,1)=2.*asin(2.*xc(i,1))+4.*xc(i,1).*sqrt(1-(4.*xc(i,1).^2));

p2xc(i,1)=sqrt(1-(4.*xc(i,1).^2));

p1x1(i,1)=2.*asin(2.*x1)+4.*x1.*sqrt(1-(4.*x1.^2));

Qstar(i,1)=(p1xc(i,1).^(3/2))./(16.*sqrt(2).*p2xc(i,1).^(1/2));

Ac(i,1)=(D.^2/8).*p1xc(i,1);

Pc(i,1)=D.*(1+asin(2.*xc(i,1)));

Rc(i,1)=Ac(i,1)./Pc(i,1);

sc(i,1)=(n.^2.*Q(i,1).^2)./(Ac(i,1).^2.*Rc(i,1).^(4/3));

f=((p1xc(i,1)./p1x1(i,1)).^(10/3).*((1+asin(2.*x1))./(1+asin(2.*xc(i,1)))).^(4/3))-s/sc;

end

end