You are using the function softmax_cross_entropy_with_logits which, according to Tensorflow's documentation, has the following specification for logits,

• logits: Per-label activations, typically a linear output. These activation energies are interpreted as unnormalized log probabilities.

Hence, you should pass the activations before the non-linearity application (in your case, softmax). You can fix it by doing the following,

def neural_network(data):
    hidden_L1 = {'weights': tf.Variable(tf.random_normal([784, neurons_L1])),
                'biases': tf.Variable(tf.random_normal([neurons_L1]))}

    hidden_L2 = {'weights': tf.Variable(tf.random_normal([neurons_L1, neurons_L2])),
                'biases': tf.Variable(tf.random_normal([neurons_L2]))}

    output_L = {'weights': tf.Variable(tf.random_normal([neurons_L2, num_of_classes])),
                'biases': tf.Variable(tf.random_normal([num_of_classes]))}

    L1 = tf.add(tf.matmul(data, hidden_L1['weights']), hidden_L1['biases']) #matrix multiplication
    L1 = tf.nn.relu(L1)

    L2 = tf.add(tf.matmul(L1, hidden_L2['weights']), hidden_L2['biases']) #matrix multiplication
    L2 = tf.nn.relu(L2)

    logits = tf.add(tf.matmul(L2, output_L['weights']), output_L['biases']) #matrix multiplication
    output = tf.nn.softmax(logits)

    return output, logits

Then, outside your function, you can retrieve the logits, and pass it to your loss function. If that doesn't work, you can also experiment with the learning rate parameter on your AdamOptimizer (see the documentation here).

You are using the function softmax_cross_entropy_with_logits which, according to Tensorflow's documentation, has the following specification for logits,

• logits: Per-label activations, typically a linear output. These activation energies are interpreted as unnormalized log probabilities.

Hence, you should pass the activations before the non-linearity application (in your case, softmax). You can fix it by doing the following,

def neural_network(data):
    hidden_L1 = {'weights': tf.Variable(tf.random_normal([784, neurons_L1])),
                'biases': tf.Variable(tf.random_normal([neurons_L1]))}

    hidden_L2 = {'weights': tf.Variable(tf.random_normal([neurons_L1, neurons_L2])),
                'biases': tf.Variable(tf.random_normal([neurons_L2]))}

    output_L = {'weights': tf.Variable(tf.random_normal([neurons_L2, num_of_classes])),
                'biases': tf.Variable(tf.random_normal([num_of_classes]))}

    L1 = tf.add(tf.matmul(data, hidden_L1['weights']), hidden_L1['biases']) #matrix multiplication
    L1 = tf.nn.relu(L1)

    L2 = tf.add(tf.matmul(L1, hidden_L2['weights']), hidden_L2['biases']) #matrix multiplication
    L2 = tf.nn.relu(L2)

    logits = tf.add(tf.matmul(L2, output_L['weights']), output_L['biases']) #matrix multiplication
    output = tf.nn.softmax(logits)

    return output, logits

Then, outside your function, you can retrieve the logits, and pass it to your loss function, as in the example bellow,

output, logits = neural_network(x)
loss = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=logits,
                                                             labels=y))

I remark that you may still be interested in the outputs tensor, for calculating your network's accuracy. If this substitution doesn't work, you should also experiment with the learning rate parameter on your AdamOptimizer (see the documentation here).