Help fitting implicit function with fit

I am trying to fit an implicit function j(E) by fitting f(j, E) to 0 for some given data I and V.

I know the true values parameters for this particular set of data and I know the function is good because when I use fsolve to get j, I get exactly my data. The true values for my parameters () are:

[j_a_true, j_c_true, alpha_anod_true, alpha_cath_true] = [2.41e-8, 0.0076533, 0.26129, 4.7057e-04, 0.75]

My problem is that using fit to get the five parameters returns -Inf for the R-square coefficient of determination and nothing close to the true values. It doesn't matter if I put in random starting values, or give the true values as starting values. R-square will always be -Inf. The difference is that when I start from the true values, SSE is reasonable.

I managed to track the problem to fit -> iGoodnessStructure. There, the SSE is very large and SST is 0. Playing around with the start values, sometimes I also get the error:

Error using fit>iFit (line 348)
Inf computed by model function, fitting cannot continue.
Try using or tightening upper and lower bounds on coefficients.

Error in fit (line 116)
[fitobj, goodness, output, convmsg] = iFit( xdatain, ydatain, fittypeobj, ...

Error in test (line 29)
[fitted, GoF] = fit([I, V], z, BV_fun,...

What am I doing wrong? Please find my code below. Any help is greatly appreciated.

% Define constants
global Faraday_const
Faraday_const = 96485.33212; % C/mol
global R_const
R_const = 8.3145; % J/K/mol
T = 298; % K
n_e = 4; % number of electrons
E_eq = 1.23; % V

V = data(:, 2); % E data
I = data(:, 1); % j data

% Define Start points, fit-function and fit curve
BV_fun = fittype( @(j_a, j_c, alpha_anod, alpha_cath, R, j_init, x) ...
j_a .* exp( (alpha_anod .* n_e .* Faraday_const) / (R_const .* T) .* (x  - j_init .* R - E_eq))-...
j_c .* exp( -(alpha_cath .* n_e .* Faraday_const) / (R_const .* T) .* (x  - j_init .* R - E_eq))-j_init,...
'independent', {'j_init', 'x'}, 'dependent', 'z');

% define array of 0s as target for f(j, E)
z = zeros(length(V), 1);

%define random sampling of starting conditions
j_a_start = random_number(0, 0.001);
j_c_start = random_number(0, 0.001);
alpha_anod_start = random_number(0, 1);
alpha_cath_start = random_number(0, 1);
R_start = random_number(0, 10);
start_values = [j_a_start, j_c_start, alpha_anod_start, alpha_cath_start, R_start];

% do the fit
[fitted, GoF] = fit([I, V], z, BV_fun,...
                   'StartPoint', start_values,...
                   'Lower', [0 0 0 0 0],...
                   'Upper', [1e-3 1e-3 1 1 10],...
                   'Algorithm', 'Trust-Region');

fits = [GoF.sse, GoF.rsquare, fitted.j_a, fitted.j_c,...
    fitted.alpha_anod, fitted.alpha_cath, fitted.R]

For running the code yourself, I am also including my data below:

data = [0.283027887344360	1.84173862266541
0.269960314035416	1.83167431879044
0.257661402225494	1.82168672609329
0.245362401008606	1.81169907379150
0.233832195401192	1.80171160030365
0.222302004694939	1.79172412681580
0.210771694779396	1.78173647451401
0.199241504073143	1.77167217063904
0.187711194157600	1.76168457794190
0.176181003451347	1.75169710445404
0.165419399738312	1.74170957136154
0.153889194130898	1.73172197866440
0.143127605319023	1.72173450517654
0.132366105914116	1.71167008209229
0.121604502201080	1.70175937938690
0.111611597239971	1.69169501590729
0.100850097835064	1.68170748281479
0.0908572077751160	1.67171989011764
0.0808643326163292	1.66173241662979
0.0716401413083076	1.65166811275482
0.0624159500002861	1.64168052005768
0.0531917512416840	1.63169298696518
0.0447362400591373	1.62170551347733
0.0370494090020657	1.61171792078018
0.0293625891208649	1.60173038768768
0.0224444400519133	1.59174279499054
0.0155263002961874	1.58175532150269
0.0101455198600888	1.57176772880554
0.00476473802700639	1.56178019571304
0.000921323895454407	1.55171577262878
-0.00138472404796630	1.54172829914093
-0.00369077292270958	1.53174082565308
-0.00522813806310296	1.52175317335129
-0.00599682098254561	1.51176569986343
-0.00599682098254561	1.50177822637558
-0.00676550390198827	1.49179057407379
-0.00599682098254561	1.48180298137665
-0.00676550390198827	1.47181550788879
-0.00676550390198827	1.46175108480454
-0.00676550390198827	1.45176361131668
-0.00676550390198827	1.44177610802651
-0.00676550390198827	1.43178848552704
-0.00676550390198827	1.42180101203919
-0.00676550390198827	1.41181350874901
-0.00676550390198827	1.40182588624954
-0.00676550390198827	1.39176161217690
-0.00676550390198827	1.38185079026222
-0.00676550390198827	1.37178648638725
-0.00676550390198827	1.36179901289940
-0.00676550390198827	1.35181139039993
-0.00676550390198827	1.34182388710976
-0.00753418589010835	1.33183629441261
-0.00753418589010835	1.32184879112244
-0.00907155219465494	1.31178439784050
-0.0106089198961854	1.30179689455032
-0.00830286927521229	1.29180930185318
-0.00753418589010835	1.28182179856300
-0.00753418589010835	1.27191109585762
-0.00676550390198827	1.26184670257568
-0.00676550390198827	1.25185919928551
-0.00676550390198827	1.24187160658836];

I = data(:, 1);
V = data(:, 2);

Answer 1

[Asterisk]

You don't need z, because you're not trying to fit a surface.

[Radu Bors]

Well I'm trying to fit f(j, E) to 0 because j(E) is an implicit function. Doesn't that make it a surface fit?