To make this question answerable, I have to assume that you have implemented back propagation correctly. That is, in terms of the underlying maths, the steps are in the right order, there is not any non-necessary repetition, and the values are calculated using correct formulae.
You can check whether the gradient values are correct by implementing gradient checking for weights, which in general you do by taking a random test case (any network architecture, random weights, inputs and target values), run your back propagation code to get your gradients, then for each weight in turn run the network forward, measure error, then adjust the weight by a very small delta, measure error again and use the difference in error to get an estimate for the gradient. If the estimate and your back propagation agree - approximately - for all weights, then your back propagation code is likely correct. If you are going to implement just one unit test for a NN implementation, this is the one to do.
To address the question directly: "Is there a more efficient way to approach backpropagation?", the answer is complicated.
The basic answer is "no". You absolutely have to do all the calculations involved in back propagation, there are no shortcuts.
However, most neural network frameworks will perform much better than your implementation, where you give an example of a network that has layers [39,39,39,1] training on 10,000 examples and taking 28 minutes to process a single epoch.
As you don't link your implementation, it is not possible to figure out precisely what the difference is. However, ignoring the possibility of bugs, the following things will probably apply:
Avoid calculating results using loops in an interpreted language. If you are looping through neuron or weight indices to multiple some values together and store back in an array in a high-level language, then this can be very inefficient. Most frameworks take advantage of fast loops and parallelisation that is possible if you frame the neural network processes (both forward and back propagation) in terms of matrix and vector manipulations. This is sometimes called vectorizaton, and most frameworks split this up so that the vector and matrix operations are black box routines implemented with strong attention to optimisation in a low-level language (like C), whilst passing the results around and deciding what to do with them can be done in higher-level language like Python or Matlab.
Mini-batches. You have to process forward and back propagation for each example seen by the network. But you don't have to apply weight updates for each example. It turns out that applying an update only once every so many examples is not only more efficient in terms of calculations, but also has benefits for convergence. Frameworks that support large amounts of parallel processing via GPU can also improve performance greatly from working in batches.
It is unlikely you will find the time to explore optimisation to the extent that the better resourced public frameworks already have. However, you can get quite far by taking advantage of vectorisation options built into your language - for instance using matrix operations in Matlab, or Numpy features in Python. This might be worthwhile so you can continue to work on and understand low-level detail in neural network models. However, at some point, if you want to work on practical problems instead of your own NN library, you will find yourself naturally switching to a third-party library because the optimisation problems and many useful additions have been solved for you already.
