Question: Curve fitting a function involving fsolve

Hello,

I'm trying to fit a double-diode equation using Maple.  I can get an analytic form of the fitting function in my two variables V and J of the form V(J), but I want J(V), so I'm using fsolve to extract this.  I'm then trying to fit the resulting function to a set of experimental data, but I'm getting an error about writing to an rtable.  I've tested the function, and it gives me a numeric value if I give it numeric inputs.  Is there a way I can fit it?  Code is below:

> fitfunc := proc (V, j1, j2, Rs1, Rs2, n1, n2, Rp1, Rp2) options operator, arrow; fsolve(V = J*Rs2+J*Rp2-j2*Rp2+n2*LambertW(j2*Rp2*exp(Rp2*(-J+j2)/vtherm)/vtherm)/vtherm+J*Rs1+J*Rp1+j1*Rp1-n1*LambertW(j1*Rp1*exp(Rp1*(J+j1)/vtherm)/vtherm)/vtherm, J) end proc;

> fitfunc(.2, 0.1e-9, 0.1e-9, 0.1e7, 0.1e7, 1.1, 1.1, 0.1e13, 0.1e13); #testing the function
                       1.000528814 10^-13

> fitdata := NonlinearFit(fitfunc, xdata, ydata, initialvalues = [0.1e-9, 0.1e-9, 0.1e7, 0.1e7, 1.1, 1.1, 0.1e11, 0.1e11]);
Error, (in Statistics:-NonlinearFit) could not store fsolve(HFloat(-0.4) = HFloat(2.0002e10)*J+HFloat(8362.310691833496)+HFloat(0.028423772659265743)*LambertW(HFloat(1.3971632824264355e9)*exp(-HFloat(2.5839793281653744e8)*J+HFloat(1.3971632824264355e9)))-HFloat(0.02842377265930558)*LambertW(HFloat(1.3971634985068152e9)*exp(HFloat(2.5839793281653744e8)*J+HFloat(1.3971634985068152e9))), J)+HFloat(8.65e-10) in a floating-point rtable