I would like to compute a Jacobian for one of my problem for a PDE, but I don’t have a function linking the unknown to my function. So I only have 2 matrices, 1 with the values of the function and 1 with the values of the unknowns.

From my shallow understanding of the problem, I cannot use automatic differentiation, because I would need the expression linking the function with the unknowns. Should I then go for finite differences to obtain a pure numerical jacobian?

I have this simple flux function, in 1D to simplify:

I’m not sure I understand your problem perfectly, but it sounds like you can use automatic differentiation.
It is not necessary to work with a “function” in the mathematical sense, that is, an explicit formula linking input to output. All you need is a “function” in the programming sense, that is, an algorithm computing the output from the input. It can use for loops, if/else statements and pretty much anything else you need.
However, it is unclear from your question if C is a function, an array or both?

Thx for your answer. I will try to clarify a bit my question with an example:

My problem is that I don’t know the relationship between the terms in the function I want to calculate the Jacobian from.

So, for example, for a function x * sin(x), let’s try to calculate \frac{\partial (x * sin(x))}{\partial x} = sin(x) + cos(x) * x

using ForwardDiff
using LinearAlgebra
using Plots
const nx = 20
const Lx = π
const Δx = Lx / (nx-1)
x = range(0, length=nx, stop= Lx)
v = zeros(nx)
v .= sin.(x)
function flux_false(x, v=v)
x .* v
end
function flux_true(x)
x .* sin.(x)
end
flux_value = flux_false(x)
Jacobian_wrong = ForwardDiff.jacobian(flux_false, x)
Jacobian_true = ForwardDiff.jacobian(flux_true, x)
deriv= sin.(x) .+ cos.(x) .* x
plot(flux_value, label="value function")
plot!(diag(Jacobian_wrong),label="jacobian wrong")
plot!(diag(Jacobian_true),label="jacobian true")
plot!(deriv, label="true derivatives", linestyle=:dash)

with the plot:

So here, I can’t use automatic differentiation because it doesn’t know the relationship between v and x, and it is just considering that v is a constant, so \frac{\partial (x * v)}{\partial x} = v = sin(x). This is the case for my real case problem.

The actual problem is that you do not include the definition of v in the flux_false function. If instead you would write

function first_step(x)
return sin.(x)
end
function second_step(x, v)
return x .* v
end
function flux(x)
v = first_step(x)
y = second_step(x, v)
return y
end

then everything would work smoothly. But your current flux_false treats v as an outside constant, with respect to which it doesn’t try to differentiate.