Something like this:

> n:= 5;
  k:= (m,l) -> (m^(1/3) + l^(1/3))^2;
  des:= {seq(diff(y[l](t),t) = 
    1/2*add(k(m,l-m)*y[m](t)*y[l-m](t), m=1..l-1)
    -y[l](t)*add(k(l,m)*y[m](t),m=1..n), l=1..n)};
  ics:= {seq(y[j](0)=1,j=1..n)};
  S:= dsolve(des union ics, numeric);
  plots[odeplot](S, [seq([t,y[j](t)],j=1..n)], t=0..2);

I don't understand: n is a constant, why would you want it to be one of the axes? Perhaps you mean e.g.

plots[odeplot](S, [seq([t,j,y[j](t)],j=1..n)], t = 0 .. 2,
  axes = box, labels = [t, j, y[j]], orientation = [250, 60]);

Dear Robert,

I hope you are doing well this sunday morning.

I am contacting you about this code I am trying to plot in 3D ( 3D surfce plot ).

Could help me please doing it ?

How can I convert this equation into matlab ?

imax:= 50; mu := 30; sigma := 5;
y[0] :=x -> exp(-(x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*Pi));
beta := (m, l) -> (m+l )^3 ;
des1:= [seq(diff(y[i](t), t) = y[i-1](t)*(sum(2^(j-i+1)*beta[i-1, j]*y[j](t), j = 1 .. i-1))+(1/2)*i*beta[i-1, i-1]*y[i-1](t)^2-y[i](t)*(sum(2^(j-i)*beta[i, j]*y[j](t), j = 1 .. i-1))-y[i](t)*(sum(beta[i, j]*y[j](t), j = 1 .. imax)), j = 1 .. imax)];
ics := {Seq(y(j)[0]=exp(-(j-mu)^2/(2*sigma^2))/(sigma*sqrt(2*Pi)); , j = 1 .. imax)};
S := dsolve(des union ics, numeric); plots[odeplot](S, [seq([t,j, y[j](t)], j = 1 ..imax)], t = 0 .. 200, axes = box, labels = [t, j, y[j]], orientation = [250, 60]);

Best regards !

sebgo