Assuming that u[0] is the same as T (if not, what is u[0]?) the following works:

eq := 4200000*(diff(T(x, t), t))-0.1e-2*(diff(T(x, t), x))^2-0.1e-2*(diff(T(x, t), x, x))*T(x, t)-.445*(diff(T(x, t), x, x))-158780+4200*T(x, t)-50000/exp(200*x);

Bc :={ T(0.03, t) = 37, D[1](T)(0,t) = 0, T(x, 0) = 38.60952381-.3512039716*cosh(64.54972242*x)-.1079529471*sinh(64.54972242*x)+.9128092390*sinh(-1.936491673+64.54972242*x) };

res:=pdsolve(eq,Bc,numeric);
res:-plot(x=0,t=0..10);
res:-animate(t=0..20,frames=5);

There are too many conditions in the first version. I'm assuming that u(t) is known. Thus from
I1(t) = I3(t) + 1/R * u(t)
follows that
I1(0) = I3(0) + 1/R* u(0)

therefore I1(0) and I3(0) cannot both be zero unless u(0) happens to be zero.

It seems that dsolve/numeric ignores the condition I1(0)=0:

restart;
daesys0:={diff(I3(t),t) = 1/L*u(t), I1(t) = I3(t) + 1/R * u(t), I3(0)=0, I1(0)=0};
u:=t->exp(t);
dsolve(daesys0,{I1(t),I3(t)}); # No solution
res:=dsolve(eval(daesys0,{R=1,L=1}),{I1(t),I3(t)},numeric);
plots:-odeplot(res,[t,I1(t)],0..1);
#Notice that I1(0) is 1 not 0.

As for your new system including the function Jv you can do as follows, taking the clue from the error message following your attempt at

res:=dsolve(daesys0,{e1(t),e2(t),Jv(t)},numeric, method =rkf45_dae,abserr=10^(-10),relerr=10^(-10),range=0..2*Pi,maxfun = 0,output=listprocedure);

restart;
with(plots):
w := 10;
u(t) :=exp(I*w*t) + 2* exp(I * 2 * w * t);
daesys0:={diff(e2(t),t) = e1(t) + Jv(t) - diff(u(t),t), 2*diff(e2(t),t) = -e2(t) - diff(u(t),t), e1(t) = -1* u(t), e1(0)=0,e2(0) = 0, Jv(0) = 50};
S:={e1=e1r+I*e1i,e2=e2r+I*e2i,Jv=Jvr+I*Jvi};
eval(daesys0,S);
s1:=(evalc@Re)~(%);
s2:=(evalc@Im)~(%%);
res:=dsolve(s1 union s2,numeric, method =rkf45_dae,abserr=10^(-10),relerr=10^(-10),range=0..2*Pi,maxfun = 0,output=listprocedure);
Xr,Yr,Zr,Xi,Yi,Zi:=op(subs(res,[e1r(t),e2r(t),Jvr(t),e1i(t),e2i(t),Jvi(t)]));
plot(Xr,0..Pi,color=red,font=[Times, roman, 14]);
plot(Xi,0..Pi,color=red,font=[Times, roman, 14]);
#odeplot doesn't need the individual procedures to be picked out:
odeplot(res,[t,e1r(t)],0..Pi);
odeplot(res,[e1r(t),e1i(t)],0..Pi,refine=1);
#If you look closer at s1 and s2 you will see that they are decoupled, i.e. the unknowns appearing in the one don't appear in the other. So you can do them separately, if you wish:
resR:=dsolve(s1,numeric, method =rkf45_dae,abserr=10^(-10),relerr=10^(-10),range=0..2*Pi,maxfun = 0,output=listprocedure);
resI:=dsolve(s2,numeric, method =rkf45_dae,abserr=10^(-10),relerr=10^(-10),range=0..2*Pi,maxfun = 0,output=listprocedure);
odeplot(resR,[t,Jvr(t)],0..Pi);

The system can be solved symbolically, but I guess you don't want that?

restart;
with(plots);
w := 10;
L :=0.001;
R :=1;
u(t) :=exp(I*w*t) + 2* exp(I * 2 * w * t);
daesys0:={diff(I3(t),t) = 1/L*u(t), I1(t) = I3(t) + 1/R * u(t), I3(0)=0};
res:=dsolve(daesys0,{I1(t),I3(t)});
X,Y:=op(subs(res,[I1(t),I3(t)]));
complexplot(X-Y,t=0..Pi/5);
complexplot(X,t=0..Pi/5,color=red);
complexplot(Y,t=0..Pi/5,color=blue);

#Notice that the condition I1(0)=0 is left out and that you need to specify the unknowns I1(t),I3(t).

If you insist on numerical work you must check your definition of dsys1. As it stands it contains equations determining t, which you surely don't want.

Try

res2:= -(diff(sol, x))/rho;

Res2 := int(rhs(res2), y = 0 .. delta(x));

#or

map(int,res2,y=0..delta(x));

#Comment: Whereas diff by itself maps over `=`, int doesn't:

diff(f(t)=g(t),t);
int(f(t)=g(t),t=a..b);
map(int,f(t)=g(t),t=a..b);

'eval' and 'subs' both need equations as in subs(a=7,expr);

Thus after having defined eq3 you can do

res:=solve({eq1, eq2, eq3}, {a, b, c});
subs(res,Te);

#Since you are evaluating Te and its derivative a few times you may find it easier to define Te as a function of y ( a procedure) and then do as follows.

restart;
Te := y->a+b*y+c*y^2;

eq1 := Te(delta[t]) = T[inf];

eq2 := D(Te)(delta[t]) = 0;

Eq3 := Te(0) = T[w]+2*(2-sigma[t])*gamma*Kn*H*D(Te)(0)/(sigma[t]*(gamma+1)*Pr);

eq3 := factor(simplify(Eq3, {2*(2-sigma[t])*gamma*Kn/(sigma[t]*(gamma+1)*Pr) = B}));

res:=solve({eq1, eq2, eq3}, {a, b, c});
subs(res,Te(y));

Below I'm using the 'known' option.

restart;
sys0 := diff(f(x), x,x)+g(x)*f(x) = 0, diff(g(x), x,x)+f(x) = 0;
bcs0 := f(0) = 1, f(1) = 2, g(0) = 0, g(1) = 1;
sol0 := dsolve({bcs0, sys0}, numeric, output = listprocedure);

#Extracting the procedures for f, f', g, g':
f0,fd0,g0,gd0:=op(subs(sol0,[f(x),diff(f(x),x),g(x),diff(g(x),x)]));

#Defining procedures for f'', g'', f''', g''':
fdd0 := proc(x) -g0(x)*f0(x) end proc:
gdd0 := proc(x) -f0(x) end proc:
fddd0 := proc(x) -g0(x)*fd0(x)-f0(x)*gd0(x) end proc:
gddd0 := proc(x) -fd0(x) end proc:

sys1 := diff(h(x),x,x)+fddd0(x)*gddd0(x) = 0, diff(k(x),x,x)-fdd0(x)-gdd0(x) = 0:
bcs1 := h(0) = 0, h(1) = 0, k(0) = 0, k(1) = 0;
#Using the known option:
sol1 := dsolve({bcs1, sys1}, numeric,known=[f0,g0,fd0,gd0,fdd0,gdd0,fddd0,gddd0],output=listprocedure);
#Extracting the procedures for h, h', k, k'
h0,hd0,k0,kd0:=op(subs(sol1,[h(x),diff(h(x),x),k(x),diff(k(x),x)]));

II haven't tried the following using components, but since 'if ... then ... end if' uses 'evalb' the following test should work on your first example.

student_answer:=(x^2-y^2)*(x^2+y^2):
if student_answer=factor(x^4-y^4) then "yes" else "no" end if;
                                    "no"

The error message from pdsolve says that an analytic solution with the given initial condition is not handled. But you can do as follows:

restart;
PDE := diff(P(x, y), y)-rho*g*beta*(T[w]-T[inf])*(1-y/delta[t](x))^2/(1+2*A*H/delta[t](x))=0;
bc := P(x, delta) = P[inf];
general := pdsolve(PDE);
subs(y = delta,bc,general);
solve(%,{_F1(x)});
sol:=subs(%,general);
#Checking the solution:
simplify(eval(PDE,sol));
simplify(eval(sol,y=delta));

How about something like this:

restart;
r:=<cos(t),sin(t),cos(2*t)>;
v:=eval(diff~(r,t),t=Pi/3);
r0:=eval(r,t=Pi/3);
T:=r0+v*s;
with(plots):
p1:=spacecurve(r,t=0..Pi/2,axes=boxed,thickness=3):
p2:=spacecurve(T,s=-1..1,axes=boxed,color=blue):
display(p1,p2);

Be aware that Maple's cylindrical coordinates are interpreted the following way (quoting the help page for plot3d,coords from Maple 12, as it happens):

"... when using cylindrical coordinates, Maple expects the command to be of the following form:
    plot3d(r(theta,z), theta=a..b, z=c..d, coords=cylindrical); "

Thus if you want the (to me) usual way z = F(r,theta), you either parametrize x,y,z by r and theta as is first done below, or you introduce z_cylindrical coordinates as is also shown below (taken from the help page again).

restart;
F:=piecewise(r=0,1,sin(6*r)/(6*r));
plot3d([r*cos(theta),r*sin(theta),F],theta=0..2*Pi,r=0..2,axes=boxed);

addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);

plot3d(F,r=0..2,theta=0..2*Pi,coords=z_cylindrical,axes=boxed);

Use isolve to get integer solutions. Then you will need to remove solutions which are not in NxN.

restart;
res:=isolve(61*a^2+1=b^2);
#The solutions as a list of lists [a,b]
RES:=map(subs,[res],[a,b]):
#The global _Z1 can take any integer value, so here are the integer solutions corresponding to _Z1 values from -3 to 3:
seq(normal(RES,expanded),_Z1=-3..3):
#However, these are not all pairs of positive integers, so:
ListTools:-FlattenOnce([%]);
select(type,%,[posint,posint]);

In any case it is best to begin by isolating the derivative, since there is some ambiguity in your case

If you want to solve the problem entirely numerically, you will have to use the parameters option as is done below.

For your differential equation it is easier to just solve the initial value problem Eq1 with y(5)=1 exactly (i.e. not numerically). And then solve the equation y(10)=5 for a (similar to what you are doing).

restart;
Eq1:=x*(diff(y(x), x))/sqrt(1+(diff(y(x), x))^2) = a;
# Eq1 implies that x*y' has the same sign as a and also that the solution only exists for x^2 > a^2.

solve(Eq1,diff(y(x),x));
#Notice below that 'isolate' only finds one solution. We can see from the conditions given that you want x>0 and y' > 0, thus a > 0. So the solution found by 'isolate' happens to be the one we want.
Eq2:=isolate(Eq1,diff(y(x),x));
dsolve({Eq2,y(5)=1});
fsolve(eval(rhs(%),x=10)=5,a);

## A purely numerical solution using parameters:
Soln := dsolve({Eq2, y(5) = 1}, numeric,parameters=[a],output=listprocedure);
Y:=subs(Soln,y(x));
p:=proc(aa) Soln(parameters=[aa]); Y(10)-5 end proc;
p(3);
A:=fsolve(p,3);
p(A);
plot(Y,A..11);

If you are satisfied with getting a finite number of coefficients, you can do like this:

eq:=(1-p)*(diff(v(p),p,p,p)+1)=H(p)*(diff(v(p),p,p,p)-25*diff(v(p),p)+1);
N:=8: #Just an example
Hs:=series(add(c[m]*p^m,m=1..N),p,N+1);
type(Hs,series);
dsolve({eval(eq,H(p)=Hs),v(0)=v0,D(v)(0)=v1,(D@D)(v)(0)=v2},v(p),type=series,order=N-2);

#Whether it should be 'order=N-2' or something else you should think about (or experiment with)
#I haven't.

If your matrix elements are functions in the mathematical sense (procedures in Maple) then the introduction of programmer indexing will cause a problem. If A is a matrix, A(2) will mean the second element, not A applied to 2.

But you can use 'apply'.

A:=Matrix([[sin],[cos]]);
apply~(A,Pi/3);
#or use map
map(apply,A,Pi/3);
#If the elements are of type 'function' in the Maple sense of 'unevaluated function call', there is no problem using subs or eval.
B:=Matrix([[sin(x)],[cos(x)]]);
subs(x=Pi/3,B);
eval(B,x=Pi/3);