Your system is linear, so we could try using the tools available for such systems.

eqs:={-4.376768352*10^5*Mea-2.542105242*10^6*Agr-1.800849667*10^6*For-60170.80926*Res+8.115163696*10^6=0,-2.542105242*10^6*Mea-1.476500135*10^7*Agr-1.045965655*10^7*For-3.494828086*10^5*Res+4.713432037*10^7=0,-1.800849667*10^6*Mea-1.045965655*10^7*Agr-7.409712508*10^6*For-2.475766891*10^5*Res+3.339036627*10^7=0,-8272.145098*Res-60170.80926*Mea-3.494828086*10^5*Agr-2.475766891*10^5*For+1.115654125*10^6=0};
A,b:=LinearAlgebra:-GenerateMatrix(eqs,[Agr,For,Mea,Res]);
LinearAlgebra:-LinearSolve(A,b);
LinearAlgebra:-GaussianElimination(A);
LinearAlgebra:-ConditionNumber(A);

You will notice that the condition number is very high. Thus you can expect numerical problems. These are seen e.g. if you assign Digits to lower values than the default (10).

The classical methods use a fixed stepsize. In the help page

?dsolve,classical

we read:
'stepsize'
Specifies the static size of each step. Note that classical does not use error correction estimates to adapt the step size. The default step size is h = min((T0-Tf)/3, 0.005) for forward integration and h = max((T0-Tf)/3, -0.005) for backward integration, respectively, where the problem is solved over the interval from T0 to Tf

#You could try inserting a procedure which remembers what it computes. Here q(t) returns 1 if t is of type numeric and t otherwise.

restart:
q:=proc(t) option remember; if not type(t,numeric) then return 'procname(_passed)' else 1 end if end proc;
#Digits:= 32;
w :=  500;
ode := diff(y(t), `$`(t, 1)) = 2*I*y(t)*q(t) + sin(w*t)*y(t)*y(t):
ics := y(0) = 1:

p2 := dsolve({ics, ode}, numeric, method = classical[rk4],output=listprocedure,known=q):
T,Y:=op(subs(p2,[t,y(t)]));
#Finding the actual step:
T(1);
R:=op(4,eval(q)):
Lix:=remove(type,map(op,[indices(R)]),name):
L:=sort(Lix);
S:=seq(L[i]-L[i-1],i=2..nops(L));
{%};

#The reason you get 0.0025 and not 0.005 is that rk4 also computes the value of f(t,y(t)) in the middle of successive steps.

Given initial conditions at zero we need to determine alpha and beta so the conditions at t=3 are satisfied.

Here is one way:
restart;
odes := diff(x(t), t) = 4*x(t)-alpha*x(t)*y(t), diff(y(t), t) = -2*y(t)+beta*x(t)*y(t);
init := dsolve({odes, x(0) = 2, y(0) = 1}, numeric, parameters = [alpha, beta], output = listprocedure);
X, Y := op(subs(init, [x(t), y(t)]));
#Check
init(parameters = [1, 1]);
plot([X, Y], 0 .. 3);
X3:=proc(a,b)
   init(parameters=[a,b]);
   eval(X)(3)-4
end proc;
Y3:=proc(a,b)
   init(parameters=[a,b]);
   eval(Y)(3)-2
end proc;
res:=fsolve({X3,Y3},[1,1]);
init(parameters=res);
X(3),Y(3);
plot([X,Y],0..3);

with(LinearAlgebra,MatrixExponential);
A:=Matrix([[1,2,0],[0,1,2],[0,0,1]]);
f:=t-><2*t*exp(t),0,0>;
x0:=<2,1,0>;
X:=MatrixExponential(A,t).x0+int~(MatrixExponential(A,t-s).f(s),s=0..t);
#Check of solution.
eval(X,t=0);
diff~(X,t) = A.X+f(t);

Your domain is B = { (s,t) | s in [0,L] and t >=0}.

The gradient of w(s,t) is the vector < diff(w(s,t),s), diff(w(s,t),t) >.
A unit normal to the part of the boundary given by (s,t) = (L,t), t>=0 is v = <1,0>.
The normal derivative is the dot product of the gradient and v, thus equal to diff(w(s,t),s).

However, your IBC-conditions contain diff(w(L,t),t) on the boundary s=L.

You can use Flatten from ListTools and then nops. If you mean that no number should be counted twice, you can use MakeUnique.

S:=[seq([seq(i, j = i .. 5)], i = 1 .. 5)];
FS:=ListTools:-Flatten(S);
nops(FS);
ListTools:-MakeUnique(FS);
nops(%);
#Obviously, if you start with a flat version you don't need Flatten:

S2:=[seq(seq(i, j = i .. 5), i = 1 .. 5)];
nops(S2);
nops(ListTools:-MakeUnique(S2));

@newlam We haven't seen your corrected version, but if you corrected the occurrence of wx(L,0)=0 to wx(L,t)=0, then it follows that D[2](wx)(L,t)=0.

There appears an H without value in IBCX and IBCY. Also there is a ws in IBCZ. Presumably these are typos.

First try the following with Digits = 10 (default), then try setting Digits to e.g. 20.

After that try an exact approach i.e. replace floating point numbers by rationals.

restart;
Digits:=20:
eqs := {2.32 = a+b, pa = .4310344828*a, pb = .4310344828*b, pa+pb = 1};
solve(eqs);
T:=LinearAlgebra:-GenerateMatrix(eqs,[a,b,pa,pb],augmented);
LinearAlgebra:-GaussianElimination(T);
LinearAlgebra:-LinearSolve(T);

restart;
eqs := convert({2.32 = a+b, pa = .4310344828*a, pb = .4310344828*b, pa+pb = 1},rational,exact);
solve(eqs);
T:=LinearAlgebra:-GenerateMatrix(eqs,[a,b,pa,pb],augmented);
LinearAlgebra:-GaussianElimination(T);
LinearAlgebra:-LinearSolve(T);

Added. The problem can be reformulated

a + b = c
s*a+s*b = 1
where s = .4310344828 and c = 2.32

This system has no solution unless 1-s*c = 0, in which case it has infinitely many.
Now s*c =1.000000000 with Digits set at 10, but 1.0000000001 with Digits set at 11, so the question is one of round-off.

In the assignment to PDEZ there appears the name wz. It ought to have been wz(s,t).

Try

indets(PDE,name);
                           {s, t, wz}
indets~([PDEX,PDEY,PDEZ],name);
                  [{s, t}, {s, t}, {s, t, wz}]

and try to run the worksheet without 'declare(....)' to find the error.

Have you tried the frames option?

The help page for odeplot:

"If the frames=n option is given, then odeplot produces an animation of the solution over the independent variable range for the plot, from lowest to highest value. The final frame of the animation is the same as the plot created without the frames option. Note: The refine option cannot be used with animations, as odeplot must choose the points for the animated plot."

You can use dchange from the PDEtools package.

The change of variables used here is  x=-2+alpha,y(x)=-1+v(alpha)

ODE := diff(y(x), x) = (y(x)+1)/(x+2)-exp((y(x)+1)/(x+2));

PDEtools:-dchange({x=-2+alpha,y(x)=-1+v(alpha)},ODE,[alpha,v]);

#You can do the whole reduction to a separable equation like this:

PDEtools:-dchange({x=-2+alpha,y(x)=-1+alpha*u(alpha)},ODE,[alpha,u]);
isolate(%,diff(u(alpha),alpha));

It seems that I ought to have seen this problem handled before, but maybe not.

IntegrationTools has an Expand that does what you expect

IntegrationTools:-Expand(Int(f(x)+g(x),x=a..b));
            
However, I don't recall seeing a similar thing for Sum (or sum).

It is easily written though. Something like this.

ES:=proc(S::specfunc(anything,{sum,Sum})) if type(op(1,S),`+`) then map(op(0,S),op(1,S),op(2,S)) end if end proc;

Then the following commands result in zero

S:=sum(sum((b[i]-b[j]),i=1..n),j=1..n);
ES(S);
simplify(%);

I don't think you will get any symbolic answer.

Replacing 'int' by 'Int' even

s1:=Int(.......)
eval(s1,{k=1,sigma=1,Q=0,theta=0});

#is not found (in Maple 15)

value(%);

#But of course you can get a numerically computed result

evalf(%);
                         0.01405865281

You are using curly braces { } where you should use square brackets [ ] so that your intended order is respected.

Try this line instead where {A,B} is replaced by [A,B].

AAA := LinearAlgebra[Eigenvectors](op(subs(s1 = 0.19e-3, [A, B]))):

It is not surprising that you get two different results when using [A,B] and [B,A].

Notice that if f is a solution to the boundary value problem, so is -f. So if there is a solution, then it is not unique. It seems that there are in fact at least 4 solutions.

You must solve the problem numerically. I didn't have any luck using dsolve/numeric directly on your problem, so I turned your problem into an initial value problem with the parameter p being D(f)(-lambda).

We must determine p s.t. f(1) = 0.

(I used the values given for beta and lambda).

init:=dsolve({s,f(-lambda)=beta*(lambda^2-1)/p,D(f)(-lambda)=p},numeric,parameters=[p],output=listprocedure);
F:=subs(init,f(x));
q:=proc(pp)
   F(parameters=[pp]);  
   F(1);
end proc;
#Testing q
q(.15);
plot(q,-.15..0.15,-1..1);
plot(q,.023..0.024,-1..1);
plot(q,.11..0.12,-1..1);
#The plots reveal approximately the values of p which make f(1) zero.

q(.0236);
q(-.0236);
F(parameters=[-.0236]);
F(1);
plot(F,-lambda..1);
F(parameters=[-.116]);
plot(F,-lambda..1);
F(1);

#An illustration

Q:=proc(pp)
   F(parameters=[pp]);  
   plot(F,-lambda..1)
end proc;
plots:-animate(Q,[p],p=-0.15..0,frames=50);