% SYSTEM OF EQUATIONS SPIN-1/2 FIELD + ELECTRIC FIELD

function S = fermions_2(t,y,mq,d,Kmin,Kmax,N,E0,w0)

j0=0;
%j0 =2*E0*w0*sech(w0.*t).^2.*tanh(w0.*t);
%j0 =-E0./(1+t).^(2);

S = zeros(N,1);

for n=0:4:N-3
    %We give values to K and W
   K=(n*d)/4 + Kmin; 
   W=(K^2+mq^2).^(1/2);
   
   %System of equations for the modes 
S(n+1) = (- (K - y(N-1)).*y(n+2) - mq*y(n+4));%Re(h1)
S(n+2) = (+ (K - y(N-1)).*y(n+1) + mq*y(n+3));%Im(h1)
S(n+3) = (+ (K - y(N-1)).*y(n+4) - mq*y(n+2));%Re(h2)
S(n+4) = (- (K - y(N-1)).*y(n+3) + mq*y(n+1));%Im(h2)

%Electric field
S(N) = S(N) - d./(2*pi)*((y(3+n).^2+y(4+n).^2) - (y(1+n).^2 + y(2+n).^2) - K./W);%pepe

end
S(N) = S(N) - y(N-1)/(pi);
S(N)=S(N) + j0;

S(N-1) = y(N);
end