I'm reading through Springer - Riemann surfaces and Farkas and Kra - Riemann surfaces and theta functions. I'm attempting to get an understanding of Weierstrass points. I've come up with a (hopefully) simple example using a hyper-elliptic function of genus 3. If I have the following curve:
$$y^2 = (x-1)(x+1)(x-2)(x+2)(x-3)(x+3)(x-4)(x+4)$$. with the following basis of holomorphic differentials.
$$\Bigg{\{}\frac{1}{y}dx, \frac{x}{y}dx, \frac{x^2}{y}dx\Bigg{\}}$$
If we were to start with the following divisor $P=(x_1,y_1)$ our basis of holomorphic differentials will consist of $\bigg{\{}\frac{(x-x_1)}{y}dx, \frac{(x-x_1)^2}{y}dx\bigg{\}}$, with dimension equal to $2$
If we were to increase $P$ to $P^2$, then the basis would change to $\bigg{\{} \frac{(x-x_1)^2}{y}dx\bigg{\}}$. Now the dimension is equal to $1$
Here's where I'm trying to make sense of things. Going from $P^2$ to $P^3$ I should not have any holomorphic differentials, but they are supposed to exists. In fact there should be $2g+2$ (8 in this case) such points. I'm thinking that in this case the Weierstrass points will be the branch points.
First of all, in the case when $X : f(x,y) = 0$ is hyperelliptic there are $2g + 2$ Weierstrass points on $X$. I believe that in general the number of points is bounded above by $g^3 - g$. Keep this in mind in case you venture into non-hyperelliptic territory.
Second, you may need to check more than just the places $P \in X$ with $x$-projections equal to one of the branch points. As a simple example, the holomorphic differentials you give are only given in affine coordinates. It may be that a place at infinity lies in the valuation divisor $(\omega_i)$ for one of the holomorphic differential basis elements with sufficient multiplicity.
The degree of the valuation divisor of a holomorphic differential $(\omega)$ is $2g-2$. Therefore, it's somewhat reasonable to expect a place $P$ occurring in $(\omega)$ with multiplicity $g$.
More concretely, it turns out that the Weierstrass points coincide with the zeros of the holomorphic $q$-differential
$$H dz^q, \quad q = g(g+1)/2$$
where
$$ H = \det \begin{pmatrix} h_1 & \cdots & h_g \\ h_1' & \cdots & h_g' \\ \vdots & & \vdots \\ h_1^{(g-1)} & \cdots & h_g^{(g-1)} \end{pmatrix} $$
and where $\omega_k = h_k(z) dz$ are local representations of a basis of holomorphic differentials. $H$ vanishes at $P$ if and only if the matrix above has a vector $(\alpha_1, \ldots, \alpha_g)^T$ in its kernel, implying that the holomorphic differential
$$\omega = \sum_{i=1}^g \alpha_i \omega_i$$
is a differential satisfying $(\omega) \geq gP$, i.e. $P$ is a Weierstrass point of $X$. (Reference: "Computational Approach to Riemann Surfaces" by Bobenko and Klein)
So I think the strategy would be to use the affine representations of the $\omega_i$'s to find the "affine Weierstrass points" and then re-parameterize at each of the places at infinity to determine if any of those are Weierstrass points as well.