Look at figure 3 (a) and (b) of the article. Each time the curve for U becomes discontinuous, a dosing is applied to the mice. If you compare this to figure 1, you will observe that the number of times dosing is applied is far greater than what figure (1) seems to suggest.
Thus what I would try if I were you is to take figure (3) and read out the time points when a jump appears and the height of the jump (which corresponds to Mi*ui in your code). This gives you two new arrays: "tdosei" and "jumpheighti". Then you can use these two arrays to repeat your computation.
If this procedure also gives results that far differ from the publication, you can additionally try to apply a simpler integration scheme with a fixed value of deltat = 0.2 days (although ode45 should of course be numerically more precise). I'm not sure whether the authors also used ode45. They report about a fixed value for the time step deltat = 0.2 days, but this could also mean that ode45 was forced to make outputs every 0.2 days (and computations were performed on a much finer grid). If the authors used a fixed time stepping scheme, this could also be a reason for results that differ from yours, and you could try such a simpler integration method.
Another thing that stroke me is that you use different arrays tDosei and tMassi. In my opinion, they must be equal because the weight of the mice was measured at the times when dosing was applied.
An indication whether you are on the right track is to check the tumor volume. It shouldn't increase far beyond 200 mm^3. In your present simulations, it reaches a maximum value of almost 1,2*10^6 mm^3 (for mouse 1) which is unrealistic.
I'd be interested if you could report about your future progress (hopefully).