restart:
with(DETools):
with(PDEtools):
u[o](r,z):=(-1/4)*diff(p[o](z),z)*(1-r^2):
ode:=gamma1*diff(u[o](r,z),z)+(1/r)*diff(v[1](r)*r,r)=0:
IC1 := v[1](0) = 0:
ans2 := dsolve({ode, IC1}, v[1](r));

The built-in command  diff  does not allow you to search for the derivative of a function by function, but if you use the well-known chain rule to differentiate the composition of functions, then we get:

q := [q1(t),q2(t),q3(t)];
L:=cos(q1(t))+sin(q2(t))+5*diff(q1(t),t) + 4*diff(q3(t),t);
diff(diff(L,t)/diff(q1(t),t), t);
                   

Edit.

h1 := solve(Vdc = 0.1500000000e-2*sqrt(2.53669508*10^8*u^3-6.06101011*10^8*u^2+3.46343435*10^8*u), u): 
A := plot([h1], Vdc = 0 .. 11.5, color = [magenta], thickness = 1): 
B := plot(Vector([0, 3.38, 5.21, 6.97, 8.4108, 10.099, 10.9232, 11.8091]), Vector([0, 0.760e-1, .1275, .1994, .2286, .3222, .3637, .999]), style = point, symbol = asterisk, color = "Blue"):
plots:-display(A, B);

For example:

restart;

eval(x+y, {x=13, y=259});

# or

assign({x=13, y=259});
x+y, x*y;
 

Should be   lambda:=3;  not  lamba:=3  in the line after  restart;

Use an inert form to avoid automatic simplification:

a%^(-3/4);
                            

restart;
plot3d(x^2 + y^2, x=0..1, y=0..x, filled=true, grid=[20,20]);

# or

plots:-shadebetween(0, x^2 + y^2, x=0..1, y=0..x, grid=[20,20]);

We see that  plots:-shadebetween   also works.
 

restart;
d1:=3: d2:=3: d3:=5:
plots:-implicitplot3d(max(x^2+y^2+z^2-d1^2, (x-3)^2+y^2+z^2-d2^2, x^2+y^2+(z-4)^2-d3^2), x=0..4, y=-3..3, z=-2..3, style=surface, color=khaki, numpoints=500000, axes=normal, scaling=constrained);
               

The user lighting was used.

Addition. We can build everything together in one plot:

restart;
d1:=3: d2:=3: d3:=5:
A:=plots:-implicitplot3d(max(x^2+y^2+z^2-d1^2,(x-3)^2+y^2+z^2-d2^2,x^2+y^2+(z-4)^2-d3^2), x=0..4,y=-3..3,z=-2..3, style=surface, color=khaki, numpoints=500000, axes=normal, scaling=constrained):    
B:=plots:-implicitplot3d([x^2+y^2+z^2=d1^2,(x-3)^2+y^2+z^2=d2^2,x^2+y^2+(z-4)^2=d3^2],x=-5..6,y=-5..5,z=-4..9,color=["LightPink","LightGreen","LightBlue"], style= wireframe, thickness=0,  numpoints=5000):
plots:-display(A, B);
                

The  roots  returns roots over the field implied by the coefficients present (see help for details).  For example, if all the coefficients are rational, then the rational roots are computed. If you want all the roots (over the field of complex numbers) use  solve  command instead of  roots :

solve(x^3-x^2-8*x+8);

Matrix(3, (i,j)->M[i, j]+b[j] );

With this method, we explicitly indicate Maple the rule for obtaining elements of the new matrix.

Another method:
M := <1,2,3; 4,5,6; 7,8,9>;
b := <10 | 11 | 12>;
M1 := <b, b, b>;
M+M1;

Replace in your code  c__1^.5  by  sqrt(c__1)  or  c__1^(1/2) . The same for  c__2^.5

 :=  instead of  =  in the first line.

restart;
g:=x->(2*x-3)^``(-3);
g(x);
g:=x->(2*x-3)^`-3`;
g(x);
 

restart;
x := (-b + sqrt(b^2 - 4*a*c) ) /(2*a):
x := subsindets(x, sqrt, t->Z);
 

In Maple, you can do animation with only one parameter. Therefore,  r  should also depend on  theta. We obtain an animation of a space curve.

An example ( a  is the parameter of the animation):

r:=sqrt(theta):
plots:-animate(plots:-spacecurve, [[r*cos(theta),r*sin(theta),r*cos(theta)], theta=0..a, color=red, thickness=2], a=0..2*Pi);

In this example, we get a spiral lying in the plane  z=x .

Edit.