As you say, the behavioral source in ngspice can only generate a voltage or current as a function of other node voltages or branch currents in your circuit. Nonetheless, you can use this, with some extra elements, to produce a nonlinear capacitor like you seem to need:
Here I used a linear CCCS and a linear capacitor to do the integration of the incoming current and track the charge variable. Then a nonlinear VCVS provides the capacitive behavior of your nonlinear capacitor. The gigohm resistor is there because SPICE requires every node to have a dc path to ground in order to obtain a solution; it will not affect the circuit solution significantly.
There's a similar example, using ngSpice syntax on page 89 in the ngSpice manual here.
Alternate solution
Also, be aware, ngSpice offers a behavioral capacitor model, described on page 71 in the manual I linked. The syntax is
CXXXXXXX n+ n C = ’expression’ <tc1 = value > < tc2 = value >
Here you should understand that the capacitance defined is the differential capacitance
\$C \equiv \frac{\mathrm{d}Q}{\mathrm{d}V}\$.
To put your equations in an appropriate form you need to re-work them a little bit. You have
\$V = K Q^2\$.
If you turn this around you have
\$Q = \sqrt{\frac{V}{K}}\$.
From which you can get
\$ \frac{\mathrm{d}Q}{\mathrm{d}V} = \frac{1}{2\sqrt{KV}}\$
Which you can easily implement in the ngSpice expression syntax.
Found the problem. Need to have "UIC" in the "tran" line.
.tran 50us 30ms 0s 500us UIC
works fine.
Guess I need to spend more time with the spice user manual.
Best Answer
Oli gave a correct answer but the I(element_name) is an extension added only to the commercial SPICE versions.
In ngspice (which is based on Berkeley Spice 3) you can only plot currents through (independent) voltage sources. These are the only currents that appear in the circuit equations SPICE works from.
In an interactive Spice session or from a special block in the script (see also this question) you can use expressions like
(v(1)-v(2))/1k
when the current is through a 1kΩ resistor between nodes 1 and 2. For reactive elements (like a 1μF capacitor) something like(v(1) - v(2))/(2*pi*frequency*1u)
should also work.