@Carl Love 

I know, but S[n] etc have huge values and for larger n the situation is worse.

Maybe the approximation is valid only near t=0. It's a kind of Picard iteration.

The "homotopy" S_n approch leads to something like:

restart;
N:=5e7;
beta:=0.00001;
mu:=0.02;
sigma:=45.6;
gamma1:=73;

S[0]:=t -> 12500000;
E[0]:=t -> 50000;
J[0]:=t -> 30000;
R[0]:=t -> 37420000;

for n from 0 to 2 do
S[n+1]:= unapply( S[n](t)-int(D(S[n])(w)- mu*N+beta*S[n](w)*J[n](w)+mu*S[n](w), w=0..t), t):
E[n+1]:= unapply( E[n](t)-int(D(E[n])(w) - beta*S[n](w)*J[n](w)+mu*E[n](w)+sigma*E[n](w), w=0..t), t);
J[n+1]:= unapply( J[n](t)-int(D(J[n])(w)- sigma*E[n](w)+mu*J[n](w)+gamma1*J[n](w), w=0..t), t);
R[n+1]:= unapply( R[n](t)-int(D(R[n])(w)-gamma1*J[n](w)+mu*R[n](w), w=0..t), t);
od;

plot(S[3], 0..70);

but does not seem to approximate the solution.

@tomleslie 

1. If you don't add the triangle inequalities, your seq test does not give a correct answer.
You can check this for  test:=[1,1,3,5,1,3];

2. The positive_definite approach is correct only for the Cayley-Menger matrix; it does not work for your  McCrea.
So the correct procedure is

fM:=proc(a,b,c,d,e,f) # a>0 etc
<0,a^2,b^2,c^2,1; a^2,0,f^2,e^2,1; b^2,f^2,0,d^2,1; c^2,e^2,d^2,0,1; 1,1,1,1,0>;
LinearAlgebra:-IsDefinite(%, query = 'positive_definite')
end;

So, it's not a matter of excessive bother.

Best regards,
V.A.

S_0(t)  etc  are needed. Probably you mean  S_0(t) = 12500000 etc.

Why don't you write your equations in Maple to be copy-pasted? Do you prefer to be typed by the forum members?

@digerdiga 

Maple 2016: For infinite sums, Maple is now more careful regarding potentially divergent parametric sums.

@digerdiga 

Mathematically the equality is not true without restrictions. The default option is now formal=false.

@nm 

But e.g. dsolve cannot solve the simple:

dsolve(diff(y(x),x)+Int((y(x)),x=0..x)=x,y(x));

and also

dsolve(diff(y(x),x)+Int((y(x)),x=0..1)=x,y(x));

@acer 

For an ODE, y(x) under Int should be (probably) also rejected.

@Adam Ledger 

1. Any programmer or designer commits mistakes ==> bugs; e,g. one of the many possible cases was not considered. Once discovered, such bugs are fixed.

2. There are theoretical aspects. E.g. there is no general algorithms able to detect whether an expression is 0 or not. If such a zero expression is not detected, a series expansion could be wrong. In this case, a workaround (of the user) could be to rewrite the expression in a simpler form.
These "bugs" cannot be fixed, but maybe a more general (not "universal" though) method could be implemented.

You should post such problems and maybe some user will come with a confirmation/workaround/explanation.

@sand15 

Any name returned from the concatenation operator will be a global variable. See ?||

@Carl Love 

I am not sure whether ad hoc is the right term. If a user asks something about a function f : R --> R,  should the answer refer to the case f : X --> Y,  with X,Y Banach spaces, in order to be non "ad hoc"?

@Carl Love 

My point is that it is not very difficult to write a procedure with a local status for a variable.
 

@Carl Love 

Maybe something like this:

restart;
x:=10:
pl:=proc(f::uneval, r::uneval)
plot(f, r);
end:
pl(sin(x), x=-Pi..Pi);

@Carl Love 

I know, my answer was in the OP's context and it is correct.