In this work, we present the numerical analysis of a bond-based peridynamic model. Model is nonlinear and is characterized by two wells. One well, near zero strain, correspond to the linear elastic interaction whereas second well corresponds to the softening of material. We present finite element analysis when the solution is in C^1 and C^2. We also analyze the finite difference approximation when the solution is in Höder space as it allows fields which do not have well-defined derivatives. We show that O(Delta t + h^2 / epsilon^2) is the rate of convergence when the solution is in C^2 and O(Delta t + h^gamma / epsilon^2) when the solution is in Höder space. h is the size of mesh, epsilon is the size of the horizon, and gamma in (0,1] is the Höder exponent. We show that the instability in the peridynamic model is due to softening of bonds and occurs when sufficiently large number of bonds have strain above critical strain. We also consider the model in one dimension. We show that the solution of peridynamic model converges to the solution of elastodynamic. The same is shown for the approximation of peridynamic model. We present some numerical results to verify the claims. This is a joint work with Dr. Robert Lipton.