I'm trying to write Physics simulation software which simulates gravitational and electric forces between n bodies.
You can have a look at my code at http://yukterez.ist.org/4kp2, sorry if I won't post the whole 1 meter long code in Latex (:
Anyway, if you have multiple forces you just add them up as vectors. For instance here the movement of particle 1 with mass m1 and charge q1 on the x axis if attracted by the particles 2 und 3 with masses m2 und m3 and the charges q2 and q3:
$a1(t) = x1''(t) =$
$$\frac{G \text{m2} (\text{x2}(t)-\text{x1}(t))}{\sqrt{\left((\text{x2}(t)-\text{x1}(t))^2+(\text{y2}(t)-\text{y1}(t))^2+(\text{z2}(t)-\text{z1}(t))^2\right)^3}}+$$
$$\frac{G \text{m3} (\text{x3}(t)-\text{x1}(t))}{\sqrt{\left((\text{x3}(t)-\text{x1}(t))^2+(\text{y3}(t)-\text{y1}(t))^2+(\text{z3}(t)-\text{z1}(t))^2\right)^3}}+$$
$$-\frac{\text{q1} \text{q2} (\text{x2}(t)-\text{x1}(t))}{((4 \pi \text{$\varepsilon $0}) \text{m1}) \sqrt{\left((\text{x2}(t)-\text{x1}(t))^2+(\text{y2}(t)-\text{y1}(t))^2+(\text{z2}(t)-\text{z1}(t))^2\right)^3}}+$$
$$-\frac{\text{q1} \text{q3} (\text{x3}(t)-\text{x1}(t))}{((4 \pi \text{$\varepsilon $0}) \text{m1}) \sqrt{\left((\text{x3}(t)-\text{x1}(t))^2+(\text{y3}(t)-\text{y1}(t))^2+(\text{z3}(t)-\text{z1}(t))^2\right)^3}}$$
Then you repeat the same procedure for y and z axis, and for the other two particles.
If you have more forces or particles just add them up.
Post Script:
But what I'm really asking is for the position as a function of time, not the acceleration.
Then just solve for the 0th derivative x(t) instead for the second one x''(t). Here is an example how you solve for acceleration, velocity and position after a given time (I choose 5 seconds, and for simplicity I reduce to 1 dimension x and 2 bodies):
Here the Objects are placed 1000 m ahead of each other on the x axis, one with mass 10000 kg and charge +1 Coulomb, the other with 5000 kg and -1 Coulomb. Initial velocity for example 10 m/s and -20 m/s, and solving for the properties of particle 1 after 5 seconds:

So the acceleration after 5 sec is 1.363 m/s², the velocity 15.45 m/s and the position 62.69 m. Understand now?
If you want to solve for particle 2 replace the x1 by x2 in the last 3 evaluations. If there are y and z axes follow the described pattern; to combine them into a plot or simulation I choose this method: Link
Post post script:
I work with Python. How do I solve it numerically?
I don't have experience with Python, but since you say
I'm struggling with the physics (and maths), not with the computer part.
I'm sure you'll manage to translate that because the Syntax is self explaining.
If you want to upgrade the math even farther for solar system simulations in this PDF, page 7 you can find the corrections for asteroid belt, solar quadrupole moment and relativity. However, you always get the terms for acceleration as second derivative which has to be solved differentially.