I think it might be a relatively trivial bug in your cost function for softmax:

J = -(sum(sum((Y).*log(h))) + lambda*p/(2*m)) 

should be

J = -sum(sum((Y).*log(h))) + lambda*p/(2*m) 

I.e. for softmax only, you have effectively subtracted the regularisation term from the cost function instead of adding it. Your back propagation calculations look correct to me if you correct this miscalculation for J.

I think it might be a relatively trivial bug in your cost function for softmax:

J = -(sum(sum((Y).*log(h))) + lambda*p/(2*m)) 

should be

J = -sum(sum((Y).*log(h)))/m + lambda*p/(2*m) 

I.e. for softmax only, you have effectively subtracted the regularisation term from the cost function instead of adding it. Also, you forgot to divide the error term by the number of examples in the batch (and you are taking this average when calculating the gradients)

Your back propagation calculations look correct to me if you correct this miscalculation for J.