It sounds like you want a parametric plot.

Example:

z:=x->ln(x+1)*(sin(x)+1);
plot([z(x),x*z(x),x=0..25],labels=[z,T]);

You could use mtaylor or D (or diff).

Example:

f:=(x,y)->sin(x)/(1+x^2)*y;
g:=(x,y)->(y-1)/(1+x^2);
f(Pi,1),g(Pi,1);
#The linalg package is deprecated so I use Jacobian from the VectorCalculus package using the long name and not doing with(VectorCalculus) as that makes a lot of confusion:
J:=VectorCalculus:-Jacobian([f(x,y),g(x,y)],[x,y]);
eval(J,{x=Pi,y=1});
mtaylor(f(x,y),[x=Pi,y=1],3);
D[1,2](f)(Pi,1);

It seems that the answer is 27/26 based on

restart;
Digits:=100:
evalf(Sum(3^(-n)*floor((1/7)*2^(n+2)), n = 1 .. infinity));
identify(%);

Proof:

Based on the familiar formula for (x^n-1)/(x-1) with x=8.

We have
S1:=floor(2^(3*k)/7)=(2^(3*k)-1)/7;
S2:=floor(2^(3*k+1)/7)=2*(2^(3*k)-1)/7;
S3:=floor(2^(3*k+2)/7)=4*(2^(3*k)-1)/7;
Sum(3^(-n)*floor((1/7)*2^(n+2)), n = 1 .. infinity)=add(Sum(3^(-3*k+2-p)*floor((1/7)*2^(3*k+p)), k = 1 .. infinity),p=0..2);
subs(S1,S2,S3,%);
value(rhs(%));

You must use a numerical approach. One that seems to work OK here is infnorm from the numapprox package.

infnorm finds the maximum of the absolute value, so you have to take care of that problem:

restart;
u:=diff(ln(cosh(t)),t$15);
Digits:=20:
M:=numapprox:-infnorm(u,t=0..50,'tm');
tm;
eval(u,t=%);
#The infnorm was attained at a value of t for which u is negative. Therefore the next step:
numapprox:-infnorm(u+M,t=0..50,'tm')-M;
tm;

A:=Matrix([[1,2],[2,3]]);
B:=<x1,x2>;
C:=<c1,c2>;
#Using elementwise equality:
eq:=A.B+C=~<0,0>;
eq[1];
eq[2];
#Alternatively:
LinearAlgebra:-GenerateEquations(A,[x1,x2],-C);

Use M[1..10,2..6] in matrixplot instead of M.

A simple example:

plot(sin,0..2*Pi);
plottools:-getdata(%);
M:=%[-1];
ExportMatrix("E:/MapleDiverse/MaplePrimes/test.dat",M);

The file test.dat can be inserted into a spreadsheet.

In the help page for DEtools[convertsys] it says:

"The initial conditions inits must each be of the form function = expression. "
A simple example:

dsolve({diff(y(x),x)=y(x), 2*y(0)=3},numeric);

This case obviously could have been handled easily also by Maple, In your case, however, convertsys would have to solve a cubic equation for (D@D)(f)(0) resulting in 3 solutions:

alpha*x+beta*x^3=f2;
solve(%,x);

Thus convertsys leaves it to the user to do the solving. Probably wise.

After with(VectorCalculus); you can do

BasisFormat(false);

The term n*(-diff(f(eta),eta$2))^(n-1)*diff(f(eta),eta$3) is likely to cause problems. Is this written exactly as you want it?

Try

evalc((-1*p)^(0.4-1)) assuming p>0;

You see that if diff(f(eta),eta,eta) is real and positive then (-diff(f(eta),eta$2))^(n-1) is imaginary.

Maple takes the principal branch of the square root.

Thus if we define s1 and s2 as

s1:=x->1/sqrt(1-x^2);
s2:=x->sqrt(1/(1-x^2));
#then there will be a difference between the following two:
s1(2);
s2(2);
#As is also seen here
1/sqrt(-3);
sqrt(-1/3);

I take it that you want a 2x182 matrix.

You can use readline:

restart;
S:="":
line := readline("E:/MapleDiverse/MaplePrimes/test.txt"):
while line <> 0 do
   S:=cat(S,line);
   line := readline("E:/MapleDiverse/MaplePrimes/test.txt")
end do:
S;
parse(S);
M:=Matrix(parse(S));

I tried the following in Maple 16:

restart;
Digits:=40:
de := proc (omega, l, M) options operator, arrow; diff(y(r), `$`(r, 2))+(2*r-2*M)*(diff(y(r), r))/(r*(r-2*M))+(omega^2*r^2/(r-2*M)^2-l*(l+1)/(r*(r-2*M)))*y(r) = 0 end proc;
init:=y(10000)=-.00008944+.00004472*I;
dinit:=D(y)(10000)=-4.4641e-6-8.9e-6*I;
sol:=dsolve([de(1/10, 1, 1), init, dinit], numeric, range = 10 .. 10000);
plots:-odeplot(sol,[Re(y(r)),Im(y(r))],10..10000,refine=1);

The plot looks fine.  Since init and dinit in your posted question are shown in an image, I didn't use all the digits.

_F5(p,q,r) is any solution to the pde mentioned inside &where, i.e.

eq2:=-_c[4]*_F5(p, q, r)+_F5(p, q, r)*lambda*p-_F5(p, q, r)*lambda+_F5(p, q, r)*gamma*r-_F5(p, q, r)*gamma+delta*(diff(_F5(p, q, r), p))-delta*(diff(_F5(p, q, r), p))*p+kappa*(diff(_F5(p, q, r), q))*N*r-kappa*(diff(_F5(p, q, r), q))*q+mu*(diff(_F5(p, q, r), r))-mu*(diff(_F5(p, q, r), r))*r+beta*(diff(_F5(p, q, r), r, p))*q*exp(-rho*tau)-beta*(diff(_F5(p, q, r), r, p))*p*r = 0;

eq2 is a pde in the variables p,q,r. Maple cannot solve eq2. You may try pdsolve(eq2);

Even if pdsolve(eq2) had produced a solution you would still have to build general solutions from the solutions obtained from the separation which took place earlier in order to satisfy initial and boundary conditions.

So you don't have much of a chance of getting any further with this problem. Numerical solution of pdes in Maple is restricted to pdes with 2 independent variables. You have 4.

You cannot use curly braces as a substitute for parentheses. Also you must use exp(x) instead of e^x.

I suppose by e^(-`ρτ`) you meant exp(-rho*tau). Thus you can try the following.

eq:=diff(G(p, q, r, t), t) = (lambda*(p-1)+gamma*(r-1))*G(p, q, r, t)+delta*(1-p)*(diff(G(p, q, r, t), p))+kappa*(N*r-q)*(diff(G(p, q, r, t), q))+mu*(1-r)*(diff(G(p, q, r, t), r))+beta*(q*exp(-rho*tau)-p*r)*(diff(G(p, q, r, t), r, p));
pdsolve(eq);

As you see, you don't get a finished solution (no surprise), but pdsolve separated variables.

You get a result that is taken a step further by adding 'build':

pdsolve(eq,build);