Fitting zeros of function to data set

What I'm trying to do here is to perform fitting of a function with 3 parameters to a data set. However, the issue at hand is the data set's dependent variable is the x-location of the minimum of the function of interest.

So, y=f(beta,x,pn) where I don't care about the value of y, but the x-location of its minimum which is a function of beta and pn. Where beta is an array of three parameters [C0,C1,R1] and pn is another parameter. Then, 0=df/dx|x=x* where x* is the location of the minimum. The anonymous function I have "dz_shift(beta,pn)" returns x* perfectly for any input "pn" (with correct beta of course). So if I have an array of pn (independent) and experimental measured x*(dependent) then I should be able to find some "beta" that can fit this data best.

What I've tried below is using nlinfit, however, it returns the error:

Error using nlinfit (line 213)

Error evaluating model function '@(beta,pn)dz_zeroed(beta(1),beta(2),beta(3),pn)/2/pi-f0'.

Caused by:

Operands to the || and && operators must be convertible to logical scalar values.

Which I'm having trouble understanding. My gut is that it is related to the fact I'm using fzero in my anonymous functions. Is there a better way to approach this problem?

clear all
close all
clc

%Initalize some constants
f0 = 4.987855028800522e+06;
c66=2.497e10;pq=2650;

%Initialize ideal data (expected)
pnm=linspace(0,1.5,16);
dfse=-f0^1.5 *sqrt(pnm/pi/c66/pq);
%experimental data (measured)
dfsm=dfse*1.3+(rand(size(pnm))-.5).*.02.*dfse;%This is x*

%Symbolic manipulation to get impedance deriviative
syms C1 R1 C0 x pn
assume(C1>0);assumeAlso(R1>0);assumeAlso(x<5.5e6*2*pi);assumeAlso(x,'real')
assumeAlso(C0>0);assumeAlso(x>4e6*2*pi);assumeAlso(pn>0);

L1=(2*pi*f0*sqrt(C1))^-2;
ws0=f0*2*pi;
L=L1+ws0*L1/pi*sqrt(2*pn/x/c66/pq);
R=R1+ws0*L1/pi*sqrt(2*pn*x/c66/pq);

%Impedance magnitude
zmag=1/(x*sqrt((C0-C1*(C1*L*x^2-1)/(C1^2*R^2*x^2+(C1*L*x^2-1)^2))^2+C1^4*R^2*x^2/(C1^2*R^2*x^2+(C1*L*x^2-1)^2)^2));
%Derivative into anon function
dz=matlabFunction(diff(zmag,x));%@(C0,C1,R1,pn,x)
%Reorganize function
dz_zeroed=@(C0,C1,R1,pnm) fzero(@(w) dz(C0,C1,R1,pnm,w),2*pi*f0);

%dz_shift returns the x-location of the minimum of zmag (x-intercept location of dz)
%ie: 0==df/dx|x=x*   where x*=dz_shift(beta,x)
dz_shift=@(beta,pn) dz_zeroed(beta(1),beta(2),beta(3),pn)/2/pi-f0;

%Nonlinear fitting
beta=[1.19e-11,26.85e-15,55.3];
%nlinfit
beta=nlinfit(pnm,dfsm,dz_shift,beta)