Solve $\frac{df(x)}{dx}=b(x)$ numerically I want to solve the following simple differential equation numerically on a grid:
$$\frac{df(x)}{dx}=b(x)$$
I tried to discretize the LHS using central differences, but the problem with that is, that the resulting linear system of equations $Ax=b$ has a singular matrix $A$ with a zero diagonal.
I also tried forward/backward difference, which gives me a (still singular) matrix $A$ with diagonal elements. I tried the Gauss-Seidel method, but got a wrong solution, because the single-sidedness of the discretization causes information to propagate only in a single direction (best seen when b is a discretized delta distribution).
My question now is:
Is there a way to discretize the equation above into a linear system $Ax=b$ which can be solved and gives the correct solution?
Thanks in advance for any help or comments.

Here's one way to use your original idea of central differences to get a numerical solution by solving a linear system.
Imagine that our domain for $f(x)$ is to be $[0,1]$, and that $f(0)=0$ is specified as the initial condition.  For a grid we will choose $n$ equally spaced points, $x_k = k/n$ for $k=0,\ldots,n$.
Now we assemble the following linear system based on the central difference approximation to the derivative:
$$ \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} = b\left( \frac{x_k + x_{k-1}}{2} \right) $$
for $k = 1,\ldots,n$.  Since $x_k - x_{k-1} = 1/n$, the above simplifies to:
$$ f(x_k) - f(x_{k-1}) = (1/n)\cdot b\left( \frac{x_k + x_{k-1}}{2} \right) $$
Now the system of unknowns $f(x_k), k=0,1,\ldots,n$ has a unique solution.  Since $f(x_0) = f(0) = 0$ is known, the remaining unknowns are expressed as summations:
$$ f(x_k) = \sum_{i=1}^k (1/n)\cdot b\left( \frac{x_i + x_{i-1}}{2} \right) $$
In other words the solution is just keeping a running sum of the indicated terms.  In a calculus course this would be called using the Midpoint Rule to approximate the integral.
Possibly you will object to evaluating the given function $b(x)$ at the midpoints of the intervals $[x_{k-1},x_k]$ on the grounds that these are not grid points.  In that case we could replace the values of $b(x)$ at the midpoints with the average of values at the endpoints of intervals:
$$ b\left( \frac{x_k + x_{k-1}}{2} \right) \approx \frac{b(x_k) + b(x_{k-1})}{2} $$
With this substitution in the system of linear equations we obtain the numerical quadrature scheme known as the Trapezoid Rule.  Since only right-hand side of the linear system is changed by this replacement, our sketch of an explicit solution is still valid with the necessary substitution inside the running summation.
