I had misunderstood your problem. When you said
I thought you meant that the dependence of (\phi, p_x, p_y, p_z) on (\psi, \theta) was defined by an equation \mathrm{lc}(\phi, p_x, p_y, p_z, \psi, \theta) = 0.
But apparently that’s not the case. Instead, you want to differentiate \mathrm{lc}(\phi(\psi, \theta), p_x(\psi, \theta), p_y(\psi, \theta), p_z(\psi, \theta), \psi, \theta) with respect to \theta, but you only have access to the functions \phi and friends through a “black box”, which you call on a grid.
What does this black box look like? Is it coded in Julia? Can you call it outside of the grid?