Ah, I'd played with this a little bit before.
Ultimately, I believe Wescott's right. Intuitively, you're emulating a resistor and diode circuit -- a system which SPICE solves iteratively, and a solution which is transcendental (you can only iterate the x[n] = e^x[n-1] form to get an approximate result, there's no closed form analytical solution -- see
formatting link
for a way of writing it).
Other sigmoid functions come to mind; you can use the Fermi-Dirac statistic (aka Logistic function), 1 / (1 + exp(-x)), to "turn on" the proportional part. But if you unwrap the tanh function, you'll see that's the same thing you're already doing (give or take some constants).
In principle, there exists a function that exactly fills in that little remainder, that you're trying to patch with the Gaussian term. It's most likely transcendental as well... and may not even have a simple description other than to being the remainder to the function in question!
Other sigmoids will have other shapes. Taking the integral of exp(-x^2) might help (erf(x)). There are other kinds out there. Do consider that, since you're SPICEing, you can do some quite excellent numerical solutions, and may not be able to express your desired relation as an equation, but you might be able to with reasonable numerical stability and few additional nodes.
Tim