Let $\frak{g}$ be a (semi-)simple lie algebra. Let $\lambda$ be a dominant integral weight. Denote $L(\lambda)$ to be the irreducible representation of highest weight $\lambda$. From BGG resolution, we know that there is a exact sequence: $$0\to M(\omega_0\cdot \lambda)\to...\to\bigoplus_{w\in W,l(w)=1} M(w\cdot \lambda)\to M(\lambda)\to L(\lambda)\to 0.$$ It is shown that $\text{dim}(\text{Hom}(M(w\cdot \lambda),M(w'\cdot \lambda)))\le 1$ and is $1$ iff $w'\le w$ in Bruhat order. When $\text{dim}(\text{Hom}(M(w\cdot \lambda),M(w'\cdot \lambda)))= 1$, we also have $M(w\cdot \lambda)$ is embedded into $M(w'\cdot \lambda)$.
As a result, we can fix a unique embedding of $M(w\cdot \lambda)$ to $M(\lambda)$. Every Verma module is generated by a highest weight vector so it makes sense to denote $v_{w\cdot \lambda}$ to be a highest weight vector of $M(w\cdot \lambda)$ inside $M(\lambda)$. It is unique up to a multiplicative constant.
Are ther any standard procedure to find all $v_{w\cdot \lambda}\in M(\lambda)$? Can they be written down as a formula?
So far, I just know that $v_{s_i\cdot \lambda}= f_i^{1+\lambda(h_i)}v_\lambda$ for each simple root $\alpha_i$.
I'm guessing that you are using the dot action: $$s_i\cdot\lambda=s_i(\lambda+\rho)-\rho.$$ Note that this is in fact an action, so $s_i\cdot(s_j\cdot\lambda)=(s_is_j)\cdot\lambda$. Now, $$v_{(s_js_i)\cdot\lambda}=v_{s_j\cdot(s_i\cdot\lambda)}=f_j^{1+(s_i\cdot\lambda)(h_j)}v_{s_i\cdot\lambda}=f_j^{1+(s_i\cdot\lambda)(h_j)}f_i^{1+\lambda(h_i)}v_\lambda.$$ Now, hopefully you can generalize to $w=s_{i_1}s_{i_2}\cdots s_{i_k}$.