What's the CFL condition of the method for your problem? This is the main thing missing from the formulation of the problem and many times at the heart of instabilities of the type you discuss. As the wikipedia site says, implicit method are usually better for maintaining a reasonable CFL condition.
The CFL condition will tell you the relationship between $\Delta t$ and $\Delta x$ that causes spurious oscillations to die down. By figuring this out in advance you can find a method that gives you a relationship that is better for your problem.
The standard way to find the CFL condition is to look for a wave solution $e^{\omega t + k\dot x}$, plug it into the numerical method and see what the numerical dispersion relation (between $\omega$ and $k$) is. The point is that for a $k$ that corresponds to an oscillation every $2\Delta x$ ($k=\pi/\Delta x$) you must have decay (meaning $Re(\omega)<0$, imaginary part not important).
The condition that the real part is negative is the CFL condition. Now, for your problem, wave solutions are probably not going to be solutions....because of the nonconstant f(x). so I would replace both by constants just to see how the sizes of f and f'' play in the CFL condition....
$\newcommand{\v}[1]{\boldsymbol{#1}}$
First about the appropriate boundary condition we would like to impose: The cross product(scalar-$\mathrm{curl}$) in $\mathbb{R}^2$ is done by we embed a vector into $\mathbb{R}^3$: $\v{v} = \langle v_1,v_2,0\rangle$, then apply the cross product, what we get would be
$$\nabla \times \v{v} = \left\langle 0,0,\left|\begin{matrix}\partial_x&\partial_y\newline v_1&v_2 \end{matrix}\right|\right\rangle
$$
Similarly the vector-$\mathbf{curl}$ is we embed a scalar function into the third dimension of $\mathbb{R}^3$: $u = \langle0,0,u\rangle$, and what we get is the rotation of the $\nabla$-operator normally denoted by $\nabla^{\perp}$:
$$
\nabla^{\perp} u = \left\langle\frac{\partial u}{\partial y},-\frac{\partial u}{\partial x},0 \right\rangle
$$
I read this way of defining 2D-curl operators in a famous computational fluid dynamics book Finite element methods for Navier-Stokes equations: theory and algorithms.
Then in $\mathbb{R}^2$ your problem can be written as:
\begin{equation}
\nabla^{\perp}(\mu^{-1}\nabla \times \v{A}) = \mu \v{J} \tag{1}
\end{equation}
Suppose we multiply a test vector $\v{v}$ vanishing on boundary on both side or the equation and integration by parts:
$$
\int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v}
- \int_{\partial \Omega} \mu^{-1}(\nabla \times \v{A})(\v{v} \times \v{n}) = \int_{\Omega} \mu\v{J}\cdot \v{v}
$$
where recall in $\mathbb{R}^2$, $B = \mu^{-1}(\nabla \times \v{A})$ is really a scalar instead of vector, hence for the tangential Neumann boundary condition, doing the $\mathbb{R}^3$-embedding we have:
$$
\v{n}\times B = \langle n_2 B, -n_1 B,0 \rangle
$$
this is the value of the scalar $B$ along the tangential direction of the domain $\v{t} = \langle n_2, -n_1,0 \rangle$, dot product with $\v{t}$ itself we could get $B$($\v{t}$ is the rotation of a unit outward normal vector), hence the Neumann type boundary condition you really wanna impose here for 2D is:
$$
\mu^{-1}(\nabla \times \v{A}) = g
$$
where $g$ is a scalar function, and your variational form becomes:
$$
\int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} = \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) \tag{2}
$$
The finite element method is we triangulating the domain, construct a finite dimensional Hilbert space $\v{V}_h$ (usually a subspace of the solution space $\v{V}$ to above variational problem), the boundary integral is treated as:
$$
\int_{\partial \Omega} g(\v{v} \cdot \v{t}) = \sum_{e\subset \partial \Omega}
\int_{e} g (\v{v}\cdot \v{t}_e)
$$
where $e$ is a boundary edge, and here is the place quadrature kicks in, that you wanna evaluate a line integral on each edge, it depends on what your finite element test function space $\v{V}_h$, nowadays a consensus is to choose Nédélec elements, in lowest order case it is $\v{v}_e = \lambda_A\nabla \lambda_B - \lambda_B\nabla \lambda_A $ for $e = \overrightarrow{AB}$.
Secondly, what is really problematic with your statement here is:
for a smooth vector field $\v{A} \in \v{C}^{\infty}(\mathbb{R}^2)$ or $\v{C}^{\infty}_c(\Omega)$, the equation (1) doesn't simply reduce to the scalar Laplace equation as you mentioned in OP, a direct observation would be that homogeneous Neumann problem for $-\Delta u = f$
has a constant as null space because a solution plus a constant is still a solution, while $\v{curl}$-problem has all the gradient vector fields as null space!
Similar thing applies for variational formulation (2), the bilinear form is not coercive, simply plugging a gradient field $\v{A} = \v{v} = \nabla \varphi$ into it. What we numerical analyst usually do here is to find a solution in the quotient space $\v{H}(\mathrm{curl})/\nabla H^1$ by adding a Lagrange multiplier to let the problem fall into mixed finite element(google MFEM if you would like to know more) framework:
$$
\begin{cases}
\int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v}
+ \int_{\Omega} \v{v}\cdot \nabla \varphi= \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) & \text{ for } \forall \v{v}\in \v{H}(\mathrm{curl})\newline
\int_{\Omega} \v{A} \cdot \nabla \psi = 0& \text{ for } \forall \psi\in H^1/\mathbb{R}
\end{cases}
$$
to make sure problem (2) is solvable using finite element, not only we have to use Lagrange multiplier to circumvent the "null space" problem, we also have to guarantee that the right side is in the range of the curl-operator, here plug $\v{v} = \nabla \varphi$ into (2):
$$
0=\int_{\Omega} \mu\v{J}\cdot \nabla \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t}) = -\int_{\Omega} \nabla \cdot(\mu\v{J}) \varphi + \int_{\partial \Omega} (\mu\v{J} \cdot \v{n}) \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t})
$$
in your case, $g=0$, also by the arbitrariness of $\varphi$, for eg we choose $\nabla \varphi \cdot \v{t} = 0$ for $\varphi\in H^1_0$, we would get so called "compatibility condition" in finite element method for Neumann problem:
$$
\nabla \cdot(\mu\v{J}) = 0\text{ in }\Omega\quad \text{ and }\quad \mu\v{J}\cdot \v{n} = 0\text{ on }\partial\Omega
$$
Lastly about the implementation of finite element method for this equation, if you have a pure Neumann problem with homogeneous boundary condition(derivative on boundary is zero), then you do not have to do anything during assembling of the stiffness matrix or right hand side integral.
Best Answer
You need to modify right-hand vector b of an equation Kx = b, where K is your stiffness matrix.
Here's how to do it, depending on which edge is on von Neumann boundary:
edge 1-2 (i.e. connecting local nodes 1 and 2),
edge 1-3,
edge 2-3,