If you take a particle approach, then you might want to investigate Particle-in-cell methods. In this particular approach, each point-charge is given a position, velocity, and charge in any dimension (1d, 2d, or 3d) and you evolve the particles using the Lorentz force law (1) and the position-derivative (2):
\begin{align}
m\frac{d\mathbf v}{dt}&=q\mathbf E+q\mathbf v\times\mathbf B \\
\frac{d\mathbf x}{dt}&=\mathbf v
\end{align}
What complicates this is that, in order to compute $\mathbf E$ and $\mathbf B$, you need to use the fact that you have a discrete collection of particles, with charge density $\rho(\mathbf x)$, and, from this, you need to determine the potential, $\phi(\mathbf x)$, and then the associated $\mathbf E$ & $\mathbf B$ fields. State-of-the-art simulations employing this are able to store information of about $10^9$ particles (note that this is on a supercomputing cluster, not on your average desktop).
If you want the fluid approach, then ACuriousMind is correct and you want to investigate Magnetohydrodynamics (MHD). These modify the Euler hydrodynamic equations to incorporate the dynamic effects of a magnetic field:
\begin{align}
\frac{\partial \rho}{\partial t}&=-\nabla\cdot\rho\mathbf v \\
\frac{\partial\rho\mathbf v}{\partial t}&=-\nabla\cdot\left(\rho\mathbf v\mathbf v+(p+\frac12B^2)\mathbb I\right) -\mathbf B\mathbf B\\
\frac{\partial e}{\partial t}&=-\nabla\cdot\left(\left[e+p+B^2+\mathbf B\mathbf B\right]\cdot\mathbf v\right)\\
\frac{\partial\mathbf B}{\partial t}&=-\nabla\cdot\left(\mathbf v\mathbf B-\mathbf B\mathbf v\right)
\end{align}
where the $\mathbf a\mathbf b$ is a dyadic, $\rho e=p/(\gamma-1)$ with $\gamma$ the adiabatic index is the equation of state, $\mathbb I$ the identity matrix, and $\rho$ the mass density, all other variables take their normal meaning. We also have the further condition that $\nabla\cdot\mathbf B=0$
In order to solve the MHD equations, you need to discretize the domain (break it up into discrete chunks and apply averaged information at "cell-centers"). The divergenceless magnetic field condition is also tough to enforce on a discrete grid, many papers have been devoted to maintaining this to machine precision! If you want to add in gravity, you need to add an extra term, $\rho\mathbf g$, to the momentum equation ($\partial_t\rho\mathbf v$ term)--this is often easier using the potential $\nabla^2\phi=\rho\mathbf g$, but also requires more work!.
If you actually do want to model this, note that you will likely be spending the next several months writing and testing the program to do whichever of the two above you've chosen.
The alternative is to go onto the github and find an MHD or PIC program that someone else spent months/years writing & testing and use that. I cannot stress enough that you will also need to spend a long time (probably on the order of weeks) studying how the code works (because it is not a black box) so that you can maybe trust the results.
