fzero coupled with ode45

I'm getting, the following comments while trying to solve ....

rspan = [0  10];
y0 =[0 0 0 0 0 0 0 0];
[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)

Error using fzero
Function values at the interval endpoints must differ in sign.

Error in solution>fun (line 24)
Opp = fzero(@(Op)fp(Op), [-10 10]);

Error in solution (line 3)
[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

function F = fun(r,y)
%% Parameters
X1 = 100;
X2 = 100;
M = 151;
S = 92;
R = 0.8/197.3;
m1 = 110;
g1 = 6.26;
m2 = 120;
g2 = 8.17;
m3 = 130;
g3 = 9.33;
Q = .5;
a = 7.33*10^(-4);
X11 = X1 - g1*(1-(a*g1*y(1)/2))*y(1);
X21 = X2 - g1*(1-(a*g1*y(1)/2))*y(1);
%% Evaluation of fp and fn
fp = @(Op)(1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2) - Op./tan(Op));
Opp = fzero(@(Op)fp(Op), [-10 10]);
 fn = @(On)(1 - (X21.*R) - sqrt(On.^2 + (X21.*R).^2) - On./tan(On));
Onn = fzero(@(On)fn(On),[-10 10]);
Ep1 = sqrt((Opp/R).^2 + X11^2);
En1 = sqrt((Onn/R).^2 + X21^2);
%% Functions
Wp = sqrt(Ep1.^2 - X11^2);
Wn = sqrt(En1.^2 - X21^2);
ap = Wp.*r;
an = Wn.*r;
Jp = besselj(0,ap);
Jp1 = besselj(1,ap);
Jn = besselj(0,an);
Jn1 = besselj(1,an);

   %% Component
   Kp = sqrt((Ep1.^2 + X11^2)./Ep1.^2) .* Jp .* r;
   Lp = sqrt((Ep1.^2 - X11^2)./Ep1.^2) .* Jp1 .* r;
   Kn = sqrt((En1.^2 + X21^2)./En1.^2) .* Jn .* r;
   Ln = sqrt((En1.^2 - X21^2)./En1.^2) .* Jn1 .* r;
%% Constant
Ap = Wp*R;
An = Wn*R;
NJp = besselj(0,Ap);
NJn = besselj(0,An);
Np = (S/(4*pi*R^3))^(1/2) .* (Ep1.*(Ep1-X11)*R).^(1/2)./(NJp.*(2*R.*Ep1.^2 - 2.*Ep1 + X11).^(1/2));
Nn = (M/(4*pi*R^3))^(1/2) .* (En1.*(En1-X21)*R).^(1/2)./(NJn.*(2*R.*En1.^2 - 2.*En1 + X21).^(1/2));  
p1 = ((Np.^2 .* (Kp.^2 - Lp.^2)) + (Nn.^2 .* (Kn.^2 - Ln.^2)))./(4*pi*r.^2);
p2 = ((Np.^2 .* (Kp.^2 + Lp.^2)) + (Nn.^2 .* (Kn.^2 + Ln.^2)))./(4*pi*r.^2);
p3 = ((Np.^2 .* (Kp.^2 + Lp.^2)) - (Nn.^2 .* (Kn.^2 + Ln.^2)))./(4*pi*r.^2);
p4 = (Np.^2 .* (Kp.^2 + Lp.^2))./(4*pi*r.^2);
%% Equations 
F(1) = y(2);
F(2) = - (2*y(2)/r) + (m1^2*y(1)) - (g1*(1- a*g1*y(1))*p1);
F(3) = y(4);
F(4) = - (2*y(4)/r) + (m2^2*y(3)) - (g2*p2);
F(5) = y(6);
F(6) = - (2*y(6)/r) + (m3^2*y(5)) - (g3*p3/2);
F(7) = y(8);
F(8) = - (2*y(8)/r) - (Q*p4);

F=[F(1);F(2);F(3);F(4);F(5);F(6);F(7);F(8)];
end

Error using fzero (line 274)

Function values at the interval endpoints must differ in sign.

Error in error_code>fun (line 26)

Opp = fzero(@(Op)fp(Op), [-10 10]);

Error in error_code>@(r,y)fun(r,y) (line 5)

[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in error_code (line 5)

[r,y] = ode45(@(r,y) fun(r,y), rspan, y0)

Warning: Function fields has the same name as a MATLAB builtin. We suggest you rename the function

to avoid a potential name conflict.

Please help me out to solve this problem, where fzero is coupled with ode45.

[Jan]

1. "Warning: Function fields has the same name as a MATLAB builtin. We suggest you rename the function to avoid a potential name conflict." - This is easy to solve.
2. You try to find the zeros of this function in the interval [-10, 10]:

I can understand why this fails.

[pritha]

I've changed what you've mentioned here. But it still not working. Still showing same error.

[Torsten]

How do you intend to handle that there are at least 4 zeros in the search interval ?

Disregarding the poles for tan(x) ...

[Jan]

@pritha: Please post the current code. "changed what you've mentioned here" does not allow to guess, what you have changed.

[pritha]

Sorry for late reply. As per your suggestion, I have changed the function field name. And for the tan part, I changed the equation to

fp = vpasolve(((Op./((1 - (Mp1.*R) - sqrt(Op.^2 + (Mp1.*R).^2)).*tan(Op))))==tan(Op), Op).

Another option i tried, I broke tan part to sin and cos part and multiplied to both of the sides and tried to calculate. For all of these above i have checked the corresponding range and put it in the code. I dont know where i'm doing the mistake again.

[Walter Roberson]

Mp1 is not defined in your code.

[pritha]

I'm sorry. It is X11. That means,

fp = vpasolve(((Op./((1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2)).*tan(Op))))==tan(Op), Op).

[Walter Roberson]

You still have a division by tan(Op) . That leads to an indefinite number of solutions that cannot be chosen between without additional information.

If you move out that tan(Op) to form

((Op./((1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2))))) == tan(Op)^2

then by examination we can see that Op == 0 is a solution except in the case that X11 == 0 or R == 0.

Answer 1

[Walter Roberson]

fp = @(Op)(1 - (X11.*R) - sqrt(Op.^2 + (X11.*R).^2) - Op./tan(Op));
Opp = fzero(@(Op)fp(Op), [-10 10]);

I was sure that I already did a write-up about this about two days ago, but I cannot find the relevant post.

You have a division by tan(Op) . That is going to be a problem if tan(Op) ever passes through zero. Which it does in the range [-10 10] at Op = -3*pi, -2*pi, -pi, 0, pi, 2*pi, and 3*pi .

fzero cannot handle functions that have a discontinuity in the search range.

When I looked at a previous post of yours (that I cannot find now), I suspected that you wanted to be working in degrees rather than radians, such as -10 degrees to +10 degrees. That would require using tand() instead of tan(). However it would still have a discontinuity at 0.

Yes it is true that limit(Op ./ tan(Op), op, 0) is 1 -- but MATLAB does not use limits in calculations unless you specifically ask for limits to be used symbolically.