one way would be to make beta (the susceptibility to infection coefficient) a function of time reflecting the time course of immunity uptake). You would then probably need numerical solutions.
that would be another approach- create two further populations infected after one vaccination and infected after two vaccinations. This would increase dimension of problem.
a further method would be to use differential-difference equations, possibly with variable lags. this would be more difficult.
ps I have not done anything in this area for many years.
You need to explain what you mean by "two separate and independent vaccination stages". Do you mean that your population of S is vaccinated once (which then converts them to an immune state) and then the previously vaccinated animals (now immune) are vaccinated a second time (which boosts their immunity and delays their return to susceptible state); or do you vaccinate the entire population of non-infected (ie. S + Immune state animals) in which case those that were fully susceptible become immune, but those previously immune receive a booster.
Once you have thought through exactly the way in which you intend to apply vaccination and the effects on the various sub-populations in your model, your can then start to formalise the equations or processes. I would think you could use a mathematical approach, or develop a mechanistic or agent based modelling approach, depending on your mathematical and / or computer programming skills.
If you are not using age but only compartments a model S -> V1 -> V2 with,for V1 and V2, different rates of return to S and different rates of contamination. The rest of the SIR scheme the same - numeric solutions not too difficult (Runge-Kutta probably the best), algebraic solutions more complicated. I would vote numeric.
I (class of) age is a concern, scheme is the same but equations more complicated and algebraic solution probably out of reach.