After the assignment to g1 I made the following changes:
type(U2y,polynom(anything, y));
#g2:=Int(U2y,x=0..X):
g2:=collect(U2y,y,s->Int(s,x=0..X)): #Int on the coefficients of the polynomial in y.
g22:=(1/delta)*g2:
g3:=U*Ux:
#g33:=Int(g3,x=0..X):
g33:=subs(x=X,U^2/2):
g4:=V*Uy:
type(g4,polynom(anything, y));
#g44:=Int(g4,x=0..X):
g44:=collect(g4,y,s->Int(s,x=0..X)):
P:=subs(x=X,g1)+g22-R*(g33+g44):
type(P,polynom(anything, y));
Eta:=subs(x=X,eta);
P1:=(1/Eta)*int(P,y=0..Eta): #int here, not Int
Po:=subs(X=0,P1):
d1:=(Po-P1):
a0:=subs(k=0.0,d1):
a1:=subs(k=0.1,d1):
a2:=subs(k=0.5,d1):
Then your plotting command gives:

MaplePrimes13-07-10i.mw

You should not use linalg. It was replaced by LinearAlgebra in Maple 6 (yes that is a long time ago).
Is this what you want:
A := <`&varphi;`[1, 1]|`&varphi;`[1, 2]|`&varphi;`[1, 3]|`&varphi;`[1, 4]>;
Notice the vertical bars.
Another syntax is:
A:=Matrix([`&varphi;`[1, 1],`&varphi;`[1, 2],`&varphi;`[1, 3],`&varphi;`[1, 4]]);

There are several equality signs that (I suppose) ought to be assignments (:=) instead.
In the assignment to PositiveEQU2, RT should be R*T.
Besides that, using R and R[f] (and more like that) is very likely to cause serious problems.
To see that you may try:
R := 8.3145; F := 96485.3365; T := 298.15; R[f] := 0;
whattype(eval(R));
eval(R);
I suggest changing the first name to Rg or similar. Using R[f] and the other indexed variables ought to be safe.
However, even after these changes, try
indets(BatSys,name);
I found the output to contain besides t and x also R[pN], R[pP].
R[pN] and R[pP] are surprising as they seem to be assigned numerical values.
However, try:
lprint(NegativePDE2);
and compare with
lprint(NegativePDE1);
The problem is with the 2d input. I never use it.
Certainly indets(BatSys,name) should only return {t, x}. So before that is the case, pdsolve/numeric won't work.

The procedure works if you change it to read:
p:=proc(pp2) local res,F0,F1,F2;
if not type([pp2],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({eq1=0,subs(p2=pp2,eq2)=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0.008,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2
end proc;

I added output=listprocedure, and removed the line res(parameters=[pp2]) as that doesn't belong here when the parameters option is not (and cannot be) used.

m:=4:
A:=LinearAlgebra:-BandMatrix([0$m,1,seq(a[i],i=1..m)], m,m+1 );
B:=Vector(m+1,i->b[i-1]);

Since K is allowed to depend on the a's, there are quite a lot of solutions:
restart;
with(LinearAlgebra):
K:=Matrix(3,symbol=b);
V:=Vector(3,symbol=a);
eq1:=a[1]+a[2]*a[3]^2*a[1]+a[2]*a[3]-a[1]^2;
eq2:=a[3]+a[1]^3+a[2]*a[3]+a[2]^2-a[1]^2;
eq3:=a[1]^3+a[2]^3+a[1]*a[2]*a[3]+a[1]*a[2]-a[3];
W:=<eq1,eq2,eq3>;
K.V=W; #Need to solve for K, i.e. the b's:
Equate(K.V,W);
res:=solve(%,{seq(seq(b[i,j],j=1..3),i=1..3)}); Using solve
K1:=subs(res,K);
K1.V;

The first example I don't understand: Where is the dependence on t? Are x and y supposed to be dependent on t?

The second equation is a partial differential equation and may also be written:

pde:=Diff(f(x,y), x,x,y,y)*x+Diff(f(x,y), x,y)*y = 0;
#The solver here is pdsolve:
pdsolve(pde);

Doesn't it follow from the boundary conditions at x = 0 and from the PDE (without solving it) that diff(u(x,t),x,x) and diff(v(x,t),x,x) are both zero at x = 0 for all t?
However, Robert Lopez gave an idea in the link given at the top. Introduce two new functions u2x and v2x:

PDE2 := {diff(u(x,t),t)=-1/20*v2x(x,t),diff(v(x,t),x,x)=v2x(x,t),
diff(v(x,t),t)=u2x(x,t),diff(u(x,t),x,x)=u2x(x,t)};
pds2 := pdsolve(PDE2, IBC, numeric, spacestep=1/50);
VAL:=pds2:-value(output=listprocedure);
VAL(0.5,3);
u2xVal,v2xVal:=op(subs(VAL,[u2x(x,t),v2x(x,t)]));
u2xVal(0,5);
v2xVal(0,1.2345);
plot3d(u2xVal,0..1,0..5);
plot3d(v2xVal,0..1,0..5);

In view of the fact that roots of polynomials are often nasty looking when they can be found symbolically at all, it may be worth your while rewriting the system as a first order linear system X' = A.X and use linear algebra. The solution is X = exp(A*t).X0, where X0 is  the initial vector X(0) = X0.
The point here is that if you replace A with its floating point version (i.e. use evalf(A) instead of A) then you are sure that eigenvectors and eigenvalues can be found so that MatrixExponential wil be able to return a result.
A similar device will not work when using dsolve since it turns floating point coefficients into rationals before doing any serious work.

restart;
eq1 := diff(x1(t),t,t) = -2*x1(t)+x2(t)-diff(x1(t),t) ;
eq2 := diff(x2(t),t,t) = -x2(t)+x1(t);
ics:= x1(0) = 1, x2(0) = 0, D(x1)(0) = 0, D(x2)(0) = 0;
#Turning the system into a first order system:
DEtools[convertsys]([eq1,eq2],[ics],[x1,x2],t,Y,YP);
subs(%[1],[YP[i]$i=1..4]);
A:=LinearAlgebra:-GenerateMatrix(%,[Y[i]$i=1..4])[1];
X:=<x1(t),x1p(t),x2(t),x2p(t)>; #The naming follows the output from convertsys
diff~(X,t) = A.X; #The system (just for show)
LinearAlgebra:-MatrixExponential(evalf(A),t).<1,0,0,0>; #Notice evalf(A) instead of A
res:=fnormal~(%); #Removing roundoff errors
plot([res[1],res[3]],t=0..30); #Plotting x1 and x2

I'm assuming that you have a couple of typos: k should have been k1 and in eq2 diff(x1(t),t,t) should have been diff(x2(t),t,t).
restart;
#k replaced by k1
eq1 := m*diff(x1(t), t, t)-f = 0; f := -k1*x1(t)+k2*(x2(t)-x1(t))-lambda*diff(x1(t), t);
#diff(x2(t), t, t) instead of diff(x1(t), t, t):
eq2 := m*diff(x2(t), t, t)-g = 0; g := -k2*(x2(t)-x1(t));
fcns := x1(t), x2(t);
ICS := {x1(0) = 1, x2(0) = 0, (D(x1))(0) = 0, (D(x2))(0) = 0};
sys := {eq1, eq2};
m := 1; k2 := 1; k1 := 1; lambda := 1;
sysdiff := sys union ICS;
#dsolve(sysdiff); #No luck: Arbitrary constants appearing
Sol_N := dsolve(sysdiff,{fcns},method=laplace);
res:=allvalues(value(Sol_N));
plot(subs(res,[x1(t),x2(t)]),t=0..30);

Maybe this small example using the parameters option can help.

restart;
eq:=diff(f(t),t,t)+f(t)=t/(t+p);
sol:=dsolve({eq,f(0)=1,D(f)(0)=0},numeric,parameters=[p],events=[[f(t)-1.2,halt],[f(t)+0.2,halt]]):
P:=proc(p) local res,q;
    sol(parameters=[p]);
    res:=op(subs(sol(15),[t,f(t),diff(f(t),t)]));
    q:=sol(eventstop);
    [p,res,q];
end proc:
P(1.4);
P(2.7);
P(10);

If the expression is not defined as a function (as in Kitonum's answer), then you can use diff and eval:

u:=x^2+2*x-3;
eval(diff(u,x),x=2);

You can still use diff and eval in the function case, but it is not as elegant as using D.
Continuing from u:

f:=unapply(u,x);
eval(diff(f(x),x),x=2);
#Compare with the simpler:
D(f)(2);

The problem is that you have a multivariate polynomial with floating point coefficients.
See the help page for factor:
?factor
For a univariate polynomial there is no problem:
u1:=1.2*x+3.4*x^2;
factor(u1);
#Now a polynomial in the variables x and a:
u:=1.2*a*x+3.4*x^2;
factor(u); #Doesn't factor
#Turning the floats into rationals:
factor(convert(u,rational,exact)); #Factors
evalf(%);
#Turning the floats into symbols instead:
evalindets(u,float,x->convert(x,symbol));
factor(%);
evalindets(%,symbol,parse);

The error message tells the story. So you are on your own.
Here is a very simple example which may or may not help you in your more complicated situation.

ode:=diff(x(t),t,t,t)=sin(x(t)); #Just one equation
dsolve({ode,x(0)=1,x(1)=0.5,x(2)=2},numeric); #Error as expected
#Now only 2 points, but with an additional requirement at 0, in this case D(x)(0) = a, where a is to be determined so that x(1) = 0.5:
p:=a->dsolve({ode,x(0)=1,x(2)=2,D(x)(0)=a},numeric,output=listprocedure);
p1:=proc(a) local X1; X1:=subs(p(a),x(t)); X1(1)-0.5 end proc;
#Just exploring the ground:
plots:-odeplot(p(-1),[t,x(t)],0..2);
#Now looking for a zero for p1, i.e. a value for a for which x(1) = 0.5 using a = -1 as an initial guess:
A:=fsolve(p1,-1);
p(A)(1);
plots:-odeplot(p(A),[t,x(t)],0..2);

If you move one line and modify it, I think you got it:
for A  from .1 by .01 to 3 do  
   for B  from .1 by .01 to 3 do
      l1:=(exp(B)-A)-1;
      l2:=(exp(B)-A/2)-1.5;
      if l1>0 and l2>0 then print(A,B); break end if
   end do
end do;

#Or if you don't want to copy from the screen:

S:=NULL:
for A  from .1 by .01 to 3 do  
   for B  from .1 by .01 to 3 do
      l1:=(exp(B)-A)-1;
      l2:=(exp(B)-A/2)-1.5;
      if l1>0 and l2>0 then S:=S,[A,B]; break end if
   end do
end do;
S; #A sequence of pairs [A,B]