MNIST in backpropagation issue accuracy not changing

I was trying to put up a script to solve the MNIST Digit Recognizer in kaggle

But I cant get the accuracy to go up. I dont know what am I missing here.

Can someone explain to me the issue here?

Can someone teach me how to add another layer in this script ?

I was mimicking this youtube video

This is it:

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

from matplotlib import pyplot as plt

data = pd.read_csv("./data/train.csv")

data.head()

data = np.array(data)
# m and n are nambre of rows and cols including labels
m, n = data.shape
print("m is", m, "and n is", n)
np.random.shuffle(data)

# transposing data to make labels in index 0
data_dev = data[0:1000].T

# labels 
Y_dev = data_dev[0]
# pixeles
X_dev = data_dev[1:n]

# here we extract training data as from 1000 to last
data_train = data[1000:m].T
Y_train = data_train[0]
X_train = data_train[1:n]


# init all weights and biases
# for the hidden and output layer
def init_params():
    # init randomly -0.5 to 0.5, shape 10 by 784
    w1 = np.random.randn(10, 784) - 0.5
    b1 = np.random.randn(10, 1) - 0.5
    w2 = np.random.randn(10, 10) - 0.5
    b2 = np.random.randn(10, 1) - 0.5
    return w1, b1, w2, b2

def ReLU(z):
    return np.maximum(0, z)

def softmax(z):
    return np.exp(z) / np.sum(np.exp(z))
   
def forward_prop(w1, b1, w2, b2, x):
    z1 = w1.dot(x) + b1
    a1 = ReLU(z1)
    z2 = w2.dot(a1) + b2
    a2 = softmax(z2)
    return z1, a1, z2, a2

# label to a number
def one_hot(y):
    one_hot_y = np.zeros((y.size, y.max() + 1))
    one_hot_y[np.arange(y.size), y] = 1
    return one_hot_y.T

def deriv_ReLU(z):
    return z > 0

def back_prop(z1, a1, z2, a2, w2, x, y):
    m = y.size
    one_hot_y = one_hot(y)
    dZ2 = a2 - one_hot_y
    dW2 = 1 / m * dZ2.dot(a1.T)
    dB2 = 1 / m * np.sum(dZ2)

    dZ1 = w2.T.dot(dZ2) * deriv_ReLU(z1)
    dW1 = 1 / m * dZ1.dot(x.T)
    dB1 = 1 / m * np.sum(dZ1)
    return dW1, dB1, dW2, dB2

def update_params(w1, b1, w2, b2, dW1, dB1, dW2, dB2, alpha):
    w1 = w1 - alpha * dW1
    b1 = b1 - alpha * dB1
    w2 = w2 - alpha * dW2
    b2 = b2 - alpha * dB2
    return w1, b1, w2, b2

def get_accuracy(a2, y):
    preds = np.argmax(a2, 0)
    return np.sum(preds == y) / y.size

def iterate(w1, b1, w2, b2, alpha, x, y, iteration):
    z1, a1, z2, a2 = forward_prop(w1, b1, w2, b2, x)
    dW1, dB1, dW2, dB2 = back_prop(z1, a1, z2, a2, w2, x, y)
    w1, b1, w2, b2 = update_params(w1, b1, w2, b2, dW1, dB1, dW2, dB2, alpha)
    if (iteration % 10 == 0):
            print("iter:", iteration)
            print("Accuracy:", get_accuracy(a2, y)) 
    return w1, b1, w2, b2

def gradient_d(x, y, iters, alpha):
    w1, b1, w2, b2 = init_params()
    for i in range(iters):
        w1, b1, w2, b2 = iterate(w1, b1, w2, b2, alpha, x, y, i)

    return w1, b1, w2, b2

w1, b1, w2, b2 = gradient_d(X_train, Y_train, 500, 0.1)