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.