want to solve a system of four odes over 'N' points

Hi, I am trying to solve ODEs associated with the number of electrons in the system (N). I want to solve the ODEs for j=1-N. I do not know where I am going wrong.

function ydot=fel4(t,y)

N=50;

ydot=[y(2);-2*y(3)*cos(y(1)+y(4));(1/N)*sum(cos(y(1)+y(4)));(-1/N)*sum(sin(y(1)+y(4)))];

Above are my functions

and my solution is:

clear all;

N=50;

tspan=[0 50];

lamdau=0.02;

rho=0.001;

lg=lamdau/(4*pi*rho);

zspan=tspan/lg;

options=odeset();

y0=[0;0;0.01;0];

[z,y]=ode45(@fel4,zspan,y0,options);

a=y(:,3);

a(N+1:end,:)=[];

z(N+1:end,:)=[];

a2=abs(a);

a3=a2.^2;

plot(z,a3)

I am not getting the right plots for this. Any help would be greatly appreciated. Thank you.

[Star Strider]

What is it supposed to look like?

Do any of the curves produced by this plot look like what you want?

y(N+1:end,:)=[];
figure(2)
plot(z, y)
grid