How to manually create a Z-table Z-tables are commonly found online. However, I am writing a precision program for this, and so I would like to find out how to calculate my own percentage values.
 A: Here is some python code to make your own z-table (in case you find it useful)
from scipy.integrate import quad
import numpy as np
import pandas as pd

def normalProbabilityDensity(x):
    constant = 1.0 / np.sqrt(2*np.pi)
    return(constant * np.exp((-x**2) / 2.0) )

standard_normal_table = pd.DataFrame(data = [],
                                 index = np.round(np.arange(0, 3.5, .1),2),
                                 columns = np.round(np.arange(0.00, .1, .01), 2))

for index in standard_normal_table.index:
    for column in standard_normal_table.columns:
        z = np.round(index + column, 2)
        value, _ = quad(normalProbabilityDensity, np.NINF, z)
        standard_normal_table.loc[index, column] = value

# Formatting to make the table look like a z-table 
standard_normal_table.index = standard_normal_table.index.astype(str)
standard_normal_table.columns = [str(column).ljust(4,'0') for column in standard_normal_table.columns]

I also derived the math explaining how to make it here. 
A: Z-tables are just values of the CDF
$$
F(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-y^{2}/2}dy
$$
at various points $x$. If you want to generate a Z-table, pick some $x$ points and evaluate that integral (using a numerical method of your choice).
As mentioned in the comments to your question, MATLAB (and other numerical software) already do this for you. Probably best not to reinvent the wheel.

Addendum: As @Ian mentioned, we should truncate the integral since it is on an unbounded domain. That is, we should look to compute instead
$$
F(x)\approx\frac{1}{\sqrt{2\pi}}\int_{w(x)}^{x}e^{-y^{2}/2}dy.
$$
The question is thus, how do we pick $w(x)$? Let's say that you are only interested in computing the integral up to $\epsilon$ accuracy. Therefore, it might be reasonable to pick $w(x)$ such that
$$
F(w(x))=\frac{1}{2}\text{erfc}(-w(x)/\sqrt{2})\leq\epsilon.
$$
A well-known bound for the erf function is
$$
\text{erfc}(x)\leq e^{-x^{2}}.
$$
Plugging this into the above,
$$
F(w(x))\leq\frac{1}{2}e^{-w(x)^{2}/2}
$$
and thus it follows that
$$
w(x)\leq-\sqrt{-2\log(2\epsilon)}.
$$

For example, if $\epsilon=10^{-6}$, $w(x)\leq-5.1230$ (approximately). This makes sense, as 5-sigma events should not make much of a contribution in computation.

This bound is very conservative, however, and you could/should come up with a tighter one, or one that takes into account relative instead of absolute error (also mentioned in comments).
A: @henry and @user104111 I will share the same answer as the thread here because I understand what you're saying. You don't want a software or tool to build a table but you need the formula & methods used to create the table from scratch and find the values in it.
So to find the values, you can proceed with more than one methods. You can use the Simpson's rule and approximate each individual value in a z score table for both negative and positive side. Alternatively, you can also use series approximation or numerical integration. As @whuber added in the other thread Mills Ratio works well out in the tails: see stats.stackexchange.com/questions/7200.
Hope this clears your doubts. Feel free to ask if you have any questions and I will elaborate my answer.
Disc: I'm affiliated with the site linked above
