How to calculate Weierstrass $p$-function by hand?

This lattice is symmetric under rotation by $\mu$, so by the homogeneity properties of $\wp$ and $\wp'$, we have $$ \wp(\mu z,\Lambda) = \mu \wp(z,\Lambda) \quad \wp'(\mu z,\Lambda) = \wp'(z,\Lambda) . $$ Applying this with the differential equation $$ \wp'^2 = 4\wp^3 - g_2 \wp - g_3 $$ forces $g_2$ to be $0$, so $ \{\wp(1/2),\wp(\mu/2),\wp((1+\mu)/2)\} = \{ \sqrt[3]{g_3/4} , \mu \sqrt[3]{g_3/4} , \mu^2 \sqrt[3]{g_3/4} \} $, since these are the roots of the right-hand side of the equation.

Since $\Lambda$ is symmetric under complex conjugation, we also have that $\wp$ is real on the real axis. This means that $\wp(1/2) =: e_1$ is real and positive (or we would have too many roots of $\wp$ in a period parallelogram by the symmetry of the lattice). This also means that $g_3=4e_1^3>0$.

This is enough information to write down an integral that gives $e_1$: since $w=\wp(z)$ satisfies $$ w'^2 = 4(w^3-e_1^3) , $$ we can invert this differential equation, using the initial condition that $w(z) \sim 1/z^2$ near $z=0$, to find that $$ z = \int_{w(z)}^{\infty} \frac{dt}{2\sqrt{t^3-e_1^3}} , $$ where the positive sign is taken for the square root; since everything here is real, the path of integration can be chosen as the real axis. (For general paths, this integral is multivalued since $w$ is doubly-periodic, but it can be made single-valued by taking appropriate branch cuts. None of this is important here.) Hence, putting $z=1/2$, we find that $$ 1 = \int_{e_1}^{\infty} \frac{dt}{\sqrt{t^3-e_1^3}} = \frac{1}{\sqrt{e_1}} \int_1^{\infty} \frac{du}{\sqrt{u^3-1}} . $$ The remaining integral is just a constant related to the $\Gamma$-function: putting $u=v^{-1/3}$ eventually gives $$ \int_1^{\infty} \frac{du}{\sqrt{u^3-1}} = \frac{\Gamma(1/3)^3}{2^{4/3}\pi} , $$ so $$ e_1 = \frac{\Gamma(1/3)^6}{2^{8/3}\pi^2} . $$

(Values at other division points can now found from this using the addition formulae, now that we know the values of $g_2,g_3$.)

This is just an addendum to Chappers' great answer. We may use the Poisson summation formula to express $\wp\left(\frac{1}{2}\right)$ in a couple of (apparently) different ways. Please allow me to use $\omega$ for denoting the third root of unity $e^{2\pi i/3}=\frac{1+i\sqrt{3}}{2}$.
We have

$$\begin{eqnarray*}\wp\left(\tfrac{1}{2}\right)&=&4+\sum_{(a,b)\neq(0,0)}\left[\frac{1}{\left(\frac{1}{2}-a-b\omega\right)^2}-\frac{1}{(a+b\omega)^2}\right]\\&=&\frac{2\pi^2}{3}+\sum_{b\geq 1}\sum_{a\in\mathbb{Z}}\left[\frac{1}{\left(\frac{1}{2}-a-b\omega\right)^2}-\frac{1}{(a+b\omega)^2}\right]\\&=&\frac{2\pi^2}{3}+\pi^2\sum_{b\geq 1}\left[\frac{1}{\cos^2\left(\frac{\pi}{2}b+\frac{\pi}{2}b i\sqrt{3}\right)}-\frac{1}{\sin^2\left(\frac{\pi}{2}b+\frac{\pi}{2}b i\sqrt{3}\right)}\right]\end{eqnarray*} $$ and by separating the contributions of even/odd values of $b$ we get $$\begin{eqnarray*}\frac{1}{\pi^2}\,\wp\left(\tfrac12\right)&=&\frac{2}{3}+\sum_{n\geq 1}\frac{1}{\cosh^2(\pi n\sqrt{3})}+\sum_{n\geq 1}\frac{1}{\sinh^2(\pi n\sqrt{3})}-\sum_{m\geq 0}\frac{1}{\sinh^2(\frac{2m+1}{2}\pi\sqrt{3})}-\sum_{m\geq 0}\frac{1}{\cosh^2(\frac{2m+1}{2}\pi\sqrt{3})}.\end{eqnarray*}$$ These series are related to special values for the complete elliptic integral of the first kind: see, for instance, this recent survey by Ce Xu. If your purpose is just a numerical evaluation of the LHS, you may just consider suitable partial sums in the RHS: the shown series are already pretty fast-convergent and they can be accelerated, too. For instance, by letting $q=e^{-2\pi\sqrt{3}}$ we have

$$ \sum_{n\geq 1}\frac{1}{4\sinh^2(\pi n\sqrt{3}) }=\sum_{n\geq 1}\frac{q^n}{(1-q^n)^2}=\sum_{n\geq 1}\sum_{k\geq 1}k q^{nk}=\sum_{N\geq 1}\sigma(N) q^N, $$ which essentially is the Eisenstein series $G_2(\tau)$ for $\tau=i\sqrt{3}$.