You sure messed up badly with the parenthesis - check them in the first and third equation!! The outer one in the third equation is unnecessary and they sure are quite wrong (and a closing one missing) in the first one.
After you correct this you have to add a fifth initial condition, something like v1(0)=???
Furthermore you have to switch the sides in your third equation as Mathcad prefers the derivatives at the LHS.
To speed up the calculation you should pre-calculate the first derivative of y(t), otherwise its calculated from anew for every iteration step. You may let the symbolic do the job by writing y1(t):=d/dt y(t) -> (symbolic evaluation) or you do the derivation yourself and define y1(t):= as the appropriate sum. Then use y1(t) instead of d/dt y(t).
But in your case y(t) should be a rectangular signal, so the derivation is zero anyway! (Apart from some single Dirac's. So you may replace d/dty(t) by zero.
Is there any reason why you are using that time consuming and inexact Fourier series instead of the pulse function which Luc provided? In your first post you showed a function where the duty cycle starts at t=0 and now you use another one where it starts at t=0.025. It would be easy to change Lucs function to do the same!
I did some changes (see picture below) but the first equation sure should read differently as with the function I came up with we get very high values at the end of the interval (10^13) which probably is not as expected.
BTW, your range definitions for t still are wrong!!
The first one (note the red dots) because the dots (..) are interpreted as unknown variable and the second one because you still use the unit s.
