You have more parameters than needed.
Let alpha-g be just alpha, a/k be just a, and beta+g be just beta.

Then you got:

sys:= (D(x))(t) = alpha*x(t)-a*x(t)^2-b*x(t)*y(t), (D(y))(t) = -beta*y(t)+c*x(t)*y(t);
#You could use the parameters option to dsolve:
res:=dsolve({sys,x(0)=1,y(0)=.1},numeric,parameters=[alpha,beta,a,b,c]);
# Now setting the parameters to the classical Volterra-Lotka model:

res(parameters=[alpha=1,beta=1,a=0,b=1,c=1]);

plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..15);
#Phase plane plot:
plots:-odeplot(res,[x(t),y(t)],0..8,view=[0..4,0..4]);

To repeat with other parameters just do similarly as in this other example:
res(parameters=[alpha=1,beta=1,a=0.1,b=1,c=1]);
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..80,size=[1000,200]);
plots:-odeplot(res,[x(t),y(t)],0..80,view=[0..4,0..4],numpoints=800);

You might also try dsolve with odeplot using that initial conditions can be changed just as when using the parameters option.
See  ?dsolve,numeric,interactive.

#Setting up a solution 'res' using the first list of initial conditions.
res:=dsolve({Eq1,Eq2,Eq3,Eq4,Eq5} union {initialset[1][]},numeric,range=0..140);
#This plot is boring because of the initial conditions we are using.
plots:-odeplot(res,[x(t),y(t)],0..140,thickness=3);

#Now defining a simple procedure pp that sets initial conditions at initialset[n] when called as pp(n).
pp:=proc(n::posint) local ini;
  ini:=[0,op(rhs~(initialset[n]))];
  res(initial=ini);
  plots:-odeplot(res,[x(t),y(t)],0..140,view=[5.2..7.2,0.2..1.6],refine=1)
end proc;
#Testing:
pp(140);
##Plotting them all:
plots:-display(seq(pp(n),n=1..180));

Solve numerically. Just add the argument 'numeric' at the end:

sol2:=dsolve(...what you got  , numeric);

It is better to give us text. That will make it possible to check it without having to type, which only (and maybe even not them) typists like.

Use plots:-odeplot to plot the solution. See help page for odeplot.

Notice that the variable(s) when given as the second argument must be given as a set if there are more than one.

See the help page for fsolve

?fsolve

where it says under parameters

"variables -  (optional) name or set(name); unknowns for which to solve"

Secondly, you are using comma as decimal separator. It should be period.
It is indeed a strange error message, but that is due to that comma-period problem.

PS: Don't use pictures, but post text or upload a worksheet. We can check what you are doing much better that way.

The left null space cannot be empty. The smallest it can be is {zerovector}.
That a basis for that space must be the empty set is correct.

So the statement: 'for every basis vector v in the basis {} we have v.ei=0 for all i'  is certainly true.

just as the statement "every unicorn is pink" is true simply because there are no unicorns.

Simpson uses the end points of the interval. Thus if the left endpoint is 0 then you have a problem.
Use a method that doesn't use the endpoints of the interval.

If you insist on using simpson then you could try replacing the function by

piecewise(x=0,1000,cos(2*x)/x^(1/3))

You will find, however, that it takes quite a number of partitions to get anywhere close to the correct result obtained by:

A:=Int(cos(2*x)/x^(1/3),x=0..1);
value(A);
evalf(%);
evalf(A);

A much better idea is to notice that the function behaves like x^(-1/3) near zero. The exact integral of that can be calculated easily. So do

f:=(cos(2*x)-1)/x^(1/3); #This function has the limit 0 at 0. Thus the piecewise below.
res1:=int(x^(-1/3),x=0..1);
res2:=evalf(ApproximateInt(piecewise(x=0,0,f), x = 0 .. 1, method = simpson,partition=8));
res1+res2;
Notice that with only partition=8 you have a reasonably good result.

You can use floor like this:

h:=t->exp(2*t);
a:=0: b:=2*Pi: p:=b-a:
plot(h(x-floor((x-a)/p)*p),x=-10..10,discont=true);

Your function in the picture is obviously not a 2*Pi periodic extension of h but looks more like a Pi periodic one.

If you do like this without assigning to n:

restart;
X := [0, 1, 2, 3];
A:=sum(sum(a[k]*X[i]^k, k = 0 .. n), i = 1 .. nops(X));
B:=sum(sum(a[k]*X[i]^k, i = 1 .. nops(X)), k = 0 .. n);
combine(A);

then you get the exact same thing in the 2 cases.

I suppose you had n:=2, so try

restart;
n:=2:
X := [0, 1, 2, 3];
sum(sum(a[k]*X[i]^k, k = 0 .. n), i = 1 .. nops(X));
#Taking it apart:
sum(a[k]*X[i]^k, k = 0 .. n); # Explained by b^0=1
sum(%,i=1..nops(X));
sum(a[k]*X[i]^k, i = 1 .. nops(X)); # Explained by 0^k=0
sum(%,k=0..n);

Maple makes the simplifications mentioned 0^k=0, b^0 = 1 for any b and any k not seen as zero.

#############
add works differently, but it only works for concrete values of n:

add(add(a[k]*X[i]^k, k = 0 .. n),i=1..nops(X));
add(add(a[k]*X[i]^k, i=1..nops(X)),k=0..n);

## sum is meant for symbolic summations (non-concrete upper or lower limits, add for the concrete cases as here.      

Use the right hand sides of

EQS:=solve({diff(eqp1,x),eqp2},{diff(u0(x),x),diff(delta0(x),x)});

as in

Fdiff:=unapply(subs(EQS,delta0(x)=F(x),u0(x)=F1(x),diff(delta0(x),x)),x,numeric);
Fdiff(.234);
F1diff:=unapply(subs(EQS,delta0(x)=F(x),u0(x)=F1(x),diff(u0(x),x)),x,numeric);

Notice now that below you will have to write Fdiff(x) and F1diff(x).

But your troubles are far from over. I have no idea if this is a viable way to go.

Incidentally, you start out by writing digits:=3; Notice that that should be spelled with a capital D. However, Digits = 3 is close to ridiculous.

res:=dsolve({Eq1,Eq2,Eq3,Eq4,ics},numeric);
plots:-odeplot(res,[t,x(t)],0..100);

See below about the new, but still undocumented(?) feature in Maple 18.

A sum can be considered an integral. That idea could be exploited.
Instead of writing a new package like IntegrationTools to handle sums one could use IntegrationTools itself:

A:=sum(a*u[k]+b*v[k],k=1..n);
subs(Int=sum,IntegrationTools:-Expand(subs(sum=Int,A)));

A number of the procedures in IntegrationTools are not done that easily, but a few are:

SIS:=proc(p::procedure,A::algebraic) subs({Int=Sum,int=sum},p(subs({Sum=Int,sum=int},A),_rest))  
  end proc;
SIS(IntegrationTools:-Expand,A);
SIS(IntegrationTools:-GetIntegrand,A);
SIS(IntegrationTools:-GetVariable,A);
SIS(IntegrationTools:-GetRange,A);
SIS(IntegrationTools:-GetParts,A);
SIS(IntegrationTools:-Flip,A); #Not correct
SIS(IntegrationTools:-Split,A,m); #This won't go but it is close.

You have { } around the integral. They cannot be used for anything but delimiters of sets. Replace them by ( ).

In fact parentheses are not needed, however.

With the model (curve) you got then you get what you got with CurveFitting:-LeastSquares.

You might have a look at Spline:

q:=CurveFitting:-Spline(y_data, x_data, k_):

plot([q,inverted_pairs],k_=min(y_data)..max(y_data),style=[line,point]);

#I corrected the assignment to inverted_pairs, "data" had disappeared):

inverted_pairs:=[seq([y_data[i], x_data[i]], i = 1 .. nops(y_data))];

The help page for alias states:

"Because aliases are resolved at the time of parsing the original input, they substitute literally, without regard to values of variables and expressions.  For example, alias(a[1]=exp(1)) followed by evalf(a[1]) will replace a[1] with exp(1) to give the result 2.718281828, but evalf(a[i]) will not substitute a[i] for its alias even when i=1.  This is because a[i] does not literally match a[1]."

that is the problem, which as discussed we can get around by not using alias at all, but that makes the code much slower as you observed.

I tested the following very simple idea on your last posted code.

1. place the alias loop on top, but it just has to be before assigning to f.

2. remove the writing, reading, and unassigning of f.

3. wrap parse(convert(...,string)) around your assignments to f[s][i] like this:

f[s][i]:=parse(convert(expand(z[s][i]-phi[s][i]),string));

This has the effect that the aliases are seen at parsing time.

The same effect is accomplished by the saving, unassigning and reading.

To see the effect in a small example compare these:

restart; ###Works
phi:=add(z[j]^2,j=1..2);
f:=z[1]-phi;
alias(z[1]=z[1](z[2]));
type(z[1],function);
hastype(f,function);
has(f,z[1]);
save f,"G:/test":
unassign('f'):
read "G:/test":

df_dz[1,2]:=diff(f,z[2]);
df_dz[1,2]:=subs(D=DD,convert(df_dz[1,2],D));

restart; ##Doesn't work
alias(z[1]=z[1](z[2]));
phi:=add(z[j]^2,j=1..2);
f:=z[1]-phi;
type(z[1],function);
hastype(f,function);
has(f,z[1]);
remove(hastype,f,function);
df_dz[1,2]:=diff(f,z[2]);
df_dz[1,2]:=subs(D=DD,convert(df_dz[1,2],D));

restart; ##Works, but notice the avoidance of symbolic indices in the definition of phi
alias(z[1]=z[1](z[2]));
phi:=z[1]^2+z[2]^2;
f:=z[1]-phi;
type(z[1],function);
hastype(f,function);
has(f,z[1]);
remove(hastype,f,function);
df_dz[1,2]:=diff(f,z[2]);
df_dz[1,2]:=subs(D=DD,convert(df_dz[1,2],D));

restart; #Works
alias(z[1]=z[1](z[2]));
phi:=add(z[j]^2,j=1..2);
f:=parse(convert(z[1]-phi,string));
type(z[1],function);
hastype(f,function);
has(f,z[1]);
remove(hastype,f,function);
df_dz[1,2]:=diff(f,z[2]);
df_dz[1,2]:=subs(D=DD,convert(df_dz[1,2],D));

You have not given r1, which appears in eqs2, a value. If you do that (I tried giving r1 the same value as r) the code should work.