Solving a differential equation numerically to plot particle path I'm trying to plot the evolution of a particle in an accretion disk by solving the equation
$$2X\frac{\partial X}{\partial\tau}=V_R(X,\tau)$$
where I have found $V_R$ numerically to be
$$V_R(X_i,\tau_n)=-\frac{3}{2X_i^2\sigma_{i,n}}\bigg[\frac{X_{i+1}\sigma_{i+1,n}-X_{i-1}\sigma_{i-1,n}}{2\Delta X}\bigg]$$
and I have all values of the surface density $\sigma_{i,j}$.
So I assume I need to integrate $V_R$ with respect to time $\tau$ in order to get $X=f(\tau)$
I'm confused when it comes to trying to do this numerically. The first initial condition is $X(0)=\sqrt{0.9}$.
Any hints or help would be greatly appreciated because I've been stuck on this problem for way too long.
 A: To solve this numerically you would start by rearranging to get
$$\frac{\partial X}{\partial \tau} = \frac{V_R(X,\tau)}{2X} $$
and then starting with $X(0)$, you find $X(\Delta\tau)$ using:
$$ X(\Delta\tau) = X(0) + \frac{\partial X}{\partial \tau}\Delta\tau $$
and so on
$$ X((n+1)\Delta\tau) = X(n\Delta\tau) + \frac{\partial X}{\partial \tau}\Delta\tau $$
It is important that your $X$ variable be continuous. In other words, don't store it as $i$ or $X_i$, but rather as $X$ (otherwise you may find your particle doesn't move anywhere!). And then when you evaluate $\partial X/\partial\tau$ you can either round $X$ to the nearest $X_i$ or you can interpolate.
A: No, you don't want to integrate $V_R$ with respect to $\tau$.  You want to solve the differential equation $2 X \frac{dX}{d\tau} = V_R(X,\tau)$ with initial condition $X(0) = \sqrt{0.9}$.  
Unless you have very special requirements, you shouldn't try to re-invent the wheel.  There is lots of good software for numerical solution of differential equations, and you can access it easily through Maple, Matlab, Mathematica etc.
