mnist with numpy only - cost does not approach 0

I'm doing course on deeplearning.ai by Andrew Ng and I wanted to try to write my own neural network to classify handwritten digits.

There's something wrong with my algorithm and I can't find the problem.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
def sigmoid(Z):
    return 1/(1 + np.exp(-Z))

def sigmoid_derivative(Z):
    return sigmoid(Z) * (1-sigmoid(Z))
train_data = pd.read_csv('train.csv')
X_train, y_train = train_data.drop(columns = ['label']).to_numpy(), train_data['label'].to_numpy().reshape(train_data.shape[0], 1)
X_train = X_train.T
y_train = y_train.T

def initialize_parameters(X, hidden_nodes):
    np.random.seed(3)
    parameters = {}
    output_layer = 10
    parameters['W1'] = np.random.randn(hidden_nodes, X.shape[0]) * np.sqrt(2/X.shape[0])
    parameters['b1'] = np.zeros((hidden_nodes, 1))
    parameters['W2'] = np.random.randn(output_layer, hidden_nodes) * np.sqrt(1/hidden_nodes)
    parameters['b2'] = np.zeros((output_layer, 1))
    return parameters

def forward_propagation(X, parameters):
    W1, b1, W2, b2 = parameters['W1'], parameters['b1'], parameters['W2'], parameters['b2']
    Z1 = np.dot(W1,X) + b1
    A1 = sigmoid(Z1)
    Z2 = np.dot(W2,A1) + b2
    A2 = sigmoid(Z2)
    forward_prop_cache = {
        'Z1' : Z1,
        'A1' : A1,
        'Z2' : Z2,
        'A2' : A2
    }
    return forward_prop_cache

def backpropagation(forward_prop_cache, parameters,X,Y):
    A1, A2, Z1, Z2 = forward_prop_cache['A1'], forward_prop_cache['A2'], forward_prop_cache['Z1'], forward_prop_cache['Z2']
    W2 = parameters['W2']
    m = X_train.shape[1]

    da2 = (A2 - Y) / (A2 * (1 - A2)) 
    dz2 = np.multiply(da2, sigmoid_derivative(Z2))
    dw2 = (1/m) * np.dot(dz2, A1.T)
    db2 = (1/m) * np.sum(dz2, keepdims=True, axis=1)
    da1 = np.dot(W2.T, dz2)
    dz1 = da1 * sigmoid_derivative(Z1)
    dw1 = (1/m) * np.dot(dz1, X.T)
    db1 = (1/m) * np.sum(dz1, keepdims=True, axis=1)
    derivatives = {
        'dw2':dw2,
        'db2':db2,
        'dw1':dw1,
        'db1':db1
    }

    return derivatives

def update_params(learning_rate, derivatives, parameters):
    dw2, db2, dw1, db1 = derivatives.values()
    W1 = parameters['W1']
    b1 = parameters['b1']
    W2 = parameters['W2']
    b2 = parameters['b2']
    W2 = W2 - learning_rate * dw2
    b2 = b2 - learning_rate * db2
    W1 = W1 - learning_rate * dw1
    b1 = b1 - learning_rate * db1
    return {'W1':W1,
            'b1':b1,
            'W2':W2,
            'b2':b2
           }

def cost_function(forward_prop_cache, Y):
    A2 = forward_prop_cache['A2']
    return -np.mean(Y * np.log(A2) + (1 - Y) * np.log(1 - A2))

def model(num_iterations, X, Y, learning_rate, hidden_nodes):
    cost = []
    parameters = initialize_parameters(X, hidden_nodes)
    
    for i in range(num_iterations):
        forward_prop_cache = forward_propagation(X, parameters)
        current_cost = cost_function(forward_prop_cache, Y)
        cost.append(current_cost)
        derivatives = backpropagation(forward_prop_cache, parameters, X, Y)
        parameters = update_params(learning_rate, derivatives, parameters)
    return cost, parameters, forward_prop_cache

num_iterations = 200
X = X_train[:,:1000]
Y = y_train[:,:1000]
learning_rate = 0.01
hidden_nodes = 200
cost, parameters, cache = model(num_iterations, X, Y, learning_rate, hidden_nodes)

Cost screenshot:

Cost plot:

As you can see there is something wrong with my algorithm. Would really appreciate some help. thank you.