You may try to isolate the time derivatives:

SYS:=solve(sys,{diff(pA(z,t),t),diff(pB(z,t),t),diff(pC(z,t),t)}):
#However, the system is then rather huge:
length(SYS);
#You don't get the error message you got before, but the computation takes longer than my patience allows for:

pdsolve(SYS,ibc,[pA,pB,pC],numeric,time=t,range=0..1);

In defining 'ode' in your corrected version there is no (implied) multiplication sign after R in ode[2].

In defining B there is a factor 2 missing.

I would never use 2D math input, especially for serious work like this.

From what you are saying it seems that you are assigning to y(t). If that is the case, don't, but assign the right hand side of the output from dsolve to a name like 'yt' or 'Y'.

RealDomain has its own 'limit' (try replacing ':' by ';' in with(RealDomain) ) and that apparently has a bug.

However, the equivalent version

limit((1-x)^(-1/x),x=0,left);

gives the correct result.

By debugging the procedure RealDomain:-limit with your original version limit((1-1/x)^(-x),x=-infinity), we see that the local variable 'ball' gets assigned the range

RealRange(Open(Float(-infinity)), Open(Float(-infinity)))

and that may be the problem as it results in the local variable 's' being assigned the value 'false' because the statement

is((1-1/x)^(-x),real) assuming x::RealRange(Open(Float(-infinity)), Open(Float(-infinity)));

results in 'false'.

[ Notice that the statement
RealRange(Open(Float(-infinity)), Open(Float(-infinity)));
evaluates to RealRange(-infinity, Open(Float(-infinity)))
thus the range is of the form [a,a), which is meaningless if a is finite. But a = -infinity, so? ]

Finally that s = false makes the procedure return 'undefined'.

How about something like

dsol1 :=dsolve({dsys1,ini1},numeric,method=lsode,var, abserr=1e-9, relerr=1e-8,output=Array([seq(k*th,k=0..N1)]));

And the matrix of values:

dsol1[2,1];

To see whether there is any justification for this method you could look at a much simpler problem.

Here I have taken a very simple problem, where the exact solution is easily found, so that the actual error can be found. If there is no connection in this simple case, there is no reason to expect a connection in a more complicated case.

restart;
Digits:=15:
sys:={diff(y(t),t) = y(t), y(0) = 1};
##The exact solution is just the exponential function exp:
dsolve(sys);
#Solving numerically:
res:=dsolve(sys,numeric, method =rkf45,abserr=10^(-10),relerr=10^(-10),output=listprocedure);
Y:=subs(res,y(t));
#Shortening your definition of S1 some (not relevant):
N := 1000:
S1 := seq([0.005 * k, Y(0.005 * k)-Y(0.005 * (k-1))], k = 2..N):
#Plotting S1 and the actual error in the same system:
plot([[S1],exp(t)-Y(t)],t=0..5,thickness=3,legend=["S1","Actual error"]);
#Plotting the actual error only:
plot(exp(t)-Y(t),t=0..5,legend="Actual error");

#You could try other simple problems. Conclusion: The answer is No!

In fact all we need to realize is that S1 would not be a sequence of zeros even if the exact solution was used for its computation, i.e. even if a perfect method was used!

I haven't checked whether your assertion is right or wrong, but you could use the StringTools package like this:

showstat(IntegralTransforms:-Tools:-SimplifyPower);
STR:=convert(eval(IntegralTransforms:-Tools:-SimplifyPower),string):
StringTools:-Search("j := evalf(j[-1][2]-1);",STR);
StringTools:-Search("k := evalf(j[1][1]+1)",STR);
StringTools:-Search("j := evalf(j[-1][2]-1); k := evalf(j[1][1]+1)",STR);

STRnew:=StringTools:-Substitute(STR,"j := evalf(j[-1][2]-1); k := evalf(j[1][1]+1)","j := evalf(t1[-1][2]-1); k := evalf(t1[1][1]+1)"):
P:=parse(STRnew):
showstat(P);

Then overwrite (if you really want to) IntegralTransforms:-Tools:-SimplifyPower. (I didn't try).

There is a syntax error. 'sys' is a set and 'values' is a sequence. Thus you need

dsolve(sys union {values}, .....);

After that you find that the implicit option is not available.

You could try a numerical solution.

indets(sys,name);
res:=dsolve(eval(sys,{seq(l[k]=1,k=1..3),seq(m[k]=1,k=1..3),g=10}) union {values},numeric,output=listprocedure);
plots:-odeplot(res,[seq([t,theta[k](t)],k=1..3)],0..0.338);

The warnings tell you the problem. The solutions (edited: their derivatives, not their values) with initial values ivp3 become infinite in finite time (in fact in less than 7 time units). Adding a range to y (and a small stepsize) may help as seen below, but notice that warnings are not errors, just warnings.

g:=y->(800-y)*y/(100+y)-450;
solve(g(y)singular(g(y),y);
sys1:=diff(y(t),t)=g(y(t));
with(DETools):

ivs3:=[y(0)=100,y(0)=250,y(0)=800]:

DEplot(sys1,y(t),t=0..10,ivs3,y=0..800,stepsize=.01);

#An alternative approach is to consider instead dt/dy = 1/(right hand side of sys1):
sys2:=diff(t(y),y)=1/g(y);
ivs2:=[t(100)=0,t(250)=0,t(800)=0];
DEplot(sys2,t(y),y=-100..800,ivs2,t=0..10,arrows=line);
p1:=%:
#Now transform the plot:
T:=plottools:-transform((y,t)->[t,y]):
T(p1);

The culprit seems to be 'radialstart', which apparently has problems with discrete data (not only list of lists of reals).

restart;
with(plots):
g:=t->3+4*cos(t);
polarplot(g(theta),theta=0..2*Pi, radialstart=-2);
N:=5:
t:=[seq(2*Pi*(i-1)/(N-1),i=1..N)]:
L:=map(s->[g(s),s],t):
polarplot(L);
plottools:-getdata(%);
polarplot(L,radialstart=-2);
plottools:-getdata(%);

#A workaround could be:

L2:=map(s->[g(s)+2,s],t):
polarplot(L2,tickmarks=[[seq(k=k-2,k=1..8)],default]);

You could use animate, e.g. like in these examples.

animate(implicitplot,[(x^2+y^2-1)^3-x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. d, grid=[50,50],gridrefine=2],d=-1.1..1.3,paraminfo=false);
animate(implicitplot,[(x^2+y^2-1)^3-d*x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. 1.3, grid=[50,50],gridrefine=2],d=-1..1,paraminfo=false);
animate(implicitplot,[(x^2+y^2-d)^3-x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. 1.3, grid=[50,50],gridrefine=2],d=0..1,paraminfo=false);
animate(implicitplot,[(x^2+y^2-d)^3-x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. 1.3, grid=[50,50],gridrefine=2,filledregions,coloring=[red,white],color=white,axes=none],d=0..1,paraminfo=false);

A text search in Maple for 'mutable' results in (among other things) 'Programming Guide, Chapter 4'.

Here is a numerical solution. The "matrix singular" situation is avoided basically by replacing l^4 by the name l4, but also by introducing l4 as a constant function:

eqn2 := diff(Y(x), x$4) = l4(x)*Y(x),diff(l4(x),x)=0;
bc := Y(0) = 0, D(Y)(0) = 1, D(D(Y))(0) = 0, Y(.5) = 0, D(Y)(.5) = 0;
sol2 := dsolve({bc, eqn2},{Y(x),l4(x)},numeric,output=listprocedure);
plots:-odeplot(sol2,[x,Y(x)]);
plots:-odeplot(sol2,[x,l4(x)^(1/4)]);
L2:=subs(sol2,l4(x)^(1/4));
L2(0);
L2(.5);

As Adri van der Meer points out there are too many conditions.

Assuming therefore that what you want is to determine L (l) so that the conditions are satisfied you can do like the following, where numerical solution of the ode is not used.

restart;
eqn := diff(Y(x), x$4) = l^4*Y(x);
bc1 := Y(0) = 0, D(Y)(0) = 1, D(D(Y))(0) = 0, Y(.5) = 0;
#bc := Y(0) = 0, (D(Y))(0) = 1, (D(D(Y)))(0) = 0, Y(.5) = 0, (D(Y))(.5) = 0;

sol := dsolve({bc1, eqn});

diff(sol,x);
eval(rhs(%),x=1/2);
#Numerator and denominator:
nm,dn:=(numer,denom)(normal(%));
plot(nm,l=0..10,-1..1);
#Here is one solution:
L:=fsolve(nm,l=7.9);
#Just making sure that the denominator is not also zero:
eval(dn,l=L);
#Plotting the solution:
eval(sol,l=L);
plot(rhs(%),x=0..1/2);
#Checking the derivative of Y at x = 1/2:
eval(diff(sol,x),{x=.5,l=L});

The system is not necessarily unsolvable as is evidenced in the following numerical solution where all constants (surely unrealistically) are set to 0.1.

cst:=indets({eqn1,eqn2,eqn3,eqn4},name) minus {T[water],T[air],T[plastic],T[Al]};
cst2:= `=`~(cst, 0.1);
fsolve(eval({eqn1,eqn2,eqn3,eqn4},cst2),{T[water],T[air],T[plastic],T[Al]});

# The output is  {T[Al] = -2.305000640, T[air] = 2.305000640, T[plastic] = -3.911849914, T[water] = 4.947857909}