First, please consult this tutorial about the LM minimization and this document on the bundle adjustment.
Along with many other derivative based descent algorithms, the Levenberg-Marquardt algorithm relies on the partial derivative matrix, a.k.a. Jacobian Matrix, which is the matrix of all first-order partial derivatives of a vector-valued function:
$$
\begin{align}
J_i &= \frac{\partial f(x_i, \beta)}{\beta} \\
\mathbf{J} &= \begin{bmatrix}
\dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_1}{\partial x_n}\\
\vdots & \ddots & \vdots\\
\dfrac{\partial f_m}{\partial x_1} & \cdots & \dfrac{\partial f_m}{\partial x_n} \end{bmatrix}
\end{align}
$$
Now how you compute Jacobian depends on what you actually want from the minimization. You could for example do a full bundle adjustment where Jacobian might have an involved expression. But, let me give an example of the simplest scheme I could think of.
Triangulation from multiple views is generally written in the form of minimization of the reprojection error. Here is the cost function for the two view problem:
$$
C(X) = d(x,\hat{x})^2 + d(x',\hat{x}')^2
$$
where $x$ and $x'$ are points in 2D views and $\hat{x}$ and $\hat{x}'$ are the projections of 3D points. In the matrix notation, the error function reads:
$$
\begin{pmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_m \\
\end{pmatrix} =
\begin{pmatrix}
\|x_1^a-P_a X^1\|+ \|x_1^b-P_b X^1\|\\
\vdots \\
\|x_m^a-P_a X^m\|+\|x_m^b-P_b X^m\| \\
\end{pmatrix}
$$
where $P_a, P_b$ are the projection matrices of cameras $a$ and $b$. Since we are searching for the 3D points only $\{\mathbf{X}_i\}$, the Jacobian is $Mx3$:
$$
\begin{align}
\mathbf{J} &= \begin{bmatrix}
\dfrac{\partial \epsilon_1}{\partial x} & \dfrac{\partial \epsilon_1}{\partial y} & \dfrac{\partial \epsilon_1}{\partial z}\\
\vdots & \ddots & \vdots\\
\dfrac{\partial \epsilon_m}{\partial x} & \dfrac{\partial \epsilon_m}{\partial y} & \dfrac{\partial \epsilon_m}{\partial z} \end{bmatrix}
\end{align}
$$
Of course this problem requires a good initial solution, which could be obtained from a linear closed form solution. The quality of the non-linear minimizer is highly dependent on the initialization and the quality of 2D points.
For the rest, and different parameterizations, check out the links.
Here are some remarks:
If you are not looking for re-inventing the wheel, then I would recommend Ceres Solver, which includes Automatic Differentiation. This saves the headache of Jacobians (and is much more accurate than numerical derivatives - numerical approximations are highly discouraged). Ceres also includes samples for bundle adjustment.
If you are also optimizing for the pose (camera orientation), then I do not recommend operating on the 12 parameter vector of rotation and translations. Rather choose a representation such as quaternions or the exponential map (twist coordinates). As this preserves the manifold-ness, you will get much higher quality solutions, without ambiguities. Some examples are given in the links I provided.
