# Problems with LsqFit on an ODE-based model

I have a 1D experimental dataset and an ODE-based model. The model is able to describe experimental data reasonably well by manual fitting of the parameter: The next step now is using `LsqFit`, at first trying to fit just one parameter:

``````julia> fit1 = curve_fit(model, tdata, ydata, [p0]);
julia> fit1.converged
true
julia> plot(fit1.resid)
`````` However it apparently doesn’t even try to fit:

``````julia> fit1.param == p0
true
``````

Now let’s change the initial value

``````julia> fit2 = curve_fit(model, tdata, ydata, [p0*1.5]);
julia> fit2.converged
true
julia> fit2.param == p0*1.5
true # !!
julia> plot(fit2.resid)

`````` The Jacobians are in both cases zero:

``````julia> all(fit1.jacobian .== 0)
true
julia> all(fit2.jacobian .== 0)
true
``````

Now I’m completely at loss Don’t use LsqFit.jl. We actively recommend against it because it’s generally not as stable as using Optimization.jl-based cost functions with a full optimizer.

Tomorrow is a project discussion with my boss, and it would be much better if we have this way of getting our material properties in addition/comparison to those used till now. Thus I may follow the advice and switch to another package later on, but today just trying to get results, quick and dirty.

I don’t need high precision, nor even precision estimation. Providing my own numeric differentiation for the Jacobian worked: Fit now converges to the same value for any reasonable initial parameters. Still would be nice to understand, why it didn’t work without it.

``````function j_m(t,p)
J = Array{Float64,2}(undef, length(t),length(p))
δrel = 1e-4
p0 = p
δpm = p0 * δrel/2
J[:,1] = (model(t, [p0+δpm]) .- model(t, [p0-δpm])) ./ (2δpm)
return J
end
``````