In  diff(f, x) the first argument  f  should be an expression not a function. Carl gave the best solution for your problem. If you still want to use  diff  command then  unapply  command helps you:

f:=x->x^2;

derivative:=unapply(diff(f(x),x), x);

derivative(2);

You can not write  

derivative:=x->diff(f(x), x);

derivative(2); 

because in this case Maple substitutes  x=2  before the derivative is calculated and an error occurs.

Here is a solution using a finite difference method. We use a simple formulas:  D(y)(x[i])=(y[i+1]-y[i-1])/2/h ,(D@@2)(y)(x[i])=(y[i+1]-2*y[i]+y[i-1])/h^2  and Simpson formula for calculation of the integral. This approach requires the boundary conditions. In the first example we specify the boundary conditions as the exact solution  y=exp(x)  has:

restart;

Eq:=x^2*diff(y(x), x, x)+50*x*diff(y(x), x)-35*y(x) = (1-exp(x+1))/(x+1)+(x^2+50*x-35)*exp(x)+Int(exp(x*t)*y(t), t = 0 .. 1):

n:=10:

h:=1/2/n:

Simpson:=L->h/3*(L[1]+2*add(L[3+2*k],k=0..n-2)+4*add(L[2+2*k],k=0..n-2)+L[-1]):

y[0]:=1: y[2*n]:=exp(1):

L:=[y[0],seq(y[i], i=1..2*n-1),y[2*n]]:

for i from 1 to 2*n-1 do

eq[i]:=eval(Eq,[diff(y(x),x)=(y[i+1]-y[i-1])/2/h, diff(y(x),x,x)=(y[i+1]-2*y[i]+y[i-1])/h^2,y(x)=y[i], Int(exp(x*t)*y(t), t = 0 .. 1)=evalf(Simpson([seq(exp(i/2/n*k/2/n)*L[k], k=1..2*n)])), x=i/2/n]);

od:

Sys:=convert(eq,set):

R:=solve(Sys);

plot([exp(x), [[0,y[0]],seq([i/2/n,rhs(R[i])],i=1..2*n-1),[1,y[2*n]]]], x=0..1, color=[blue,red], thickness=0, view=[0..1,0..y[2*n]]);  # Visual comparison of exact and approximate solutions (the plots merge)

In the second example we specify the other boundary conditions and get a different solution:

restart:

Eq:=x^2*diff(y(x), x, x)+50*x*diff(y(x), x)-35*y(x) = (1-exp(x+1))/(x+1)+(x^2+50*x-35)*exp(x)+Int(exp(x*t)*y(t), t = 0 .. 1):

n:=10:

h:=1/2/n:

Simpson:=L->h/3*(L[1]+2*add(L[3+2*k],k=0..n-2)+4*add(L[2+2*k],k=0..n-2)+L[-1]):

y[0]:=1: y[2*n]:=1.5:

L:=[y[0],seq(y[i], i=1..2*n-1),y[2*n]];

for i from 1 to 2*n-1 do

eq[i]:=eval(Eq,[diff(y(x),x)=(y[i+1]-y[i-1])/2/h, diff(y(x),x,x)=(y[i+1]-2*y[i]+y[i-1])/h^2,y(x)=y[i], Int(exp(x*t)*y(t), t = 0 .. 1)=evalf(Simpson([seq(exp(i/2/n*k/2/n)*L[k], k=1..2*n)])), x=i/2/n]);

od:

Sys:=convert(eq,set):

R:=solve(Sys);

plot([[0,y[0]],seq([i/2/n,rhs(R[i])],i=1..2*n-1),[1,y[2*n]]],x=0..1,color=red, thickness=0, view=[0..1,0..y[2*n]]);

Can we conclude that there is more than one solution of the original problem?

See help on  Explore  command. To plot vectors you can use  plots[arrow]  command. The imaginary unit is coded (by default) in Maple as  I  not  j .

Of course the second calculation  evalf(gg(2))=-0.3333333333  is incorrect. Probably Maple uses formal geometric series summation  b+b*q+b*q^2+..+b*q^n+..=b/(1-q)  which is correct only if  |q|<1 .

Indeed for the series  gg(2)=Sum(2^jj*2^jj, jj = 0 .. infinity)   we have  b=1 ,  q=4    and   1/(1-4)=-0.3333333333

You can search for a numerical solution using  numeric  option. To do this, all the parameters and initial conditions must be specified. With the resulting procedure, you can find solutions in some points, build plots, etc.

sys := {diff(x(t), t) = -k*x(t)*sqrt(x(t)^2+y(t)^2)/m, diff(y(t), t) = -k*y(t)*sqrt(x(t)^2+y(t)^2)/m - g}:

sol := dsolve({eval(op(sys), {g = 9.8, k = 1, m = 2}), x(0) = 0, y(0) = 0}, numeric);

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 2, color = [red, blue], thickness = 3);

If you mean the 50 roots of this equation, then here are the first 50 positive roots:

seq(fsolve(0.5*lambda*tan(lambda)-1, lambda=`if`(n=0,0..Pi/2,-Pi/2+n*Pi..Pi/2+n*Pi)), n=0..49);

 Addition. You can easily plot these points with the following code (for clarity, plotted the first 10 roots in the form of small red circles):

L:=[seq(fsolve(0.5*lambda*tan(lambda)-1, lambda=`if`(n=0,0..Pi/2,-Pi/2+n*Pi..Pi/2+n*Pi)), n=0..49)]:

plots[display](plot(0.5*lambda*tan(lambda)-1, lambda=0..Pi/2+9*Pi+0.01, -10..10, color=blue, discont), plot(map(t->[t,0], L[1..10]), style=point, symbol=solidcircle, symbolsize=12, color=red, scaling=constrained));

Complex function of a complex argument may be regarded as a mapping from R^2 into R^2. On the left - the space of arguments,  on the right - the space of their images:

R:=evalc(1/(1-(x+I*y)));

R1:=op(1,R),coeff(R,I);

f:=unapply([R1],x,y);

F:=plottools[transform](f):

A:=plot([seq([convert(<cos(Pi/4*i),-sin(Pi/4*i);sin(Pi/4*i),cos(Pi/4*i)>.<t,0>,list)[], t=-0.8..0.8], i=0..7), seq([r*cos(t),r*sin(t),t=0..2*Pi], r=0.2..0.8,0.2)], color=[red,gold,brown,yellow,blue,green,violet,cyan], thickness=3):

plots[display](<A | F(A)>, scaling=constrained);

 Edited.

If we go to polar coordinates, Maple can easily find the solution in explicit form, even simpler than Mathematica does:

restart;

T:={x(t)=r(t)*cos(phi(t)), y(t)=r(t)*sin(phi(t))}:  # The transformation to polar coordinates

sys:={diff(x(t),t) = -k/m*x(t)*sqrt(x(t)^2+y(t)^2), diff(y(t),t) = -k/m*y(t)*sqrt(x(t)^2+y(t)^2)}:  # The original system

subs(csgn(r(t))=1, simplify(eval(sys,T)));  # The transformation of original system to polar coordinates

sys1:=solve(%,{diff(r(t),t),diff(phi(t),t)});  # The original system in polar coordinates in standard form

dsolve(sys1, {r(t),phi(t)});  # The general solution in polar coordinates

dsolve(sys1 union {r(0)=1,phi(0)=Pi/3},{r(t),phi(t)});  # The solutin with some initial conditions

eval(T, %);  # The same solution in Cartesian coordinates

In fact, the point moves along a straight line toward the origin.

ConvIntoIneq  procedure solves the problem:

ConvIntoIneq:=proc(Rel)

local t, R, x, y;

t:=indets(Rel)[];

R:=solve(Rel);

x, y:=op(1,R), op(2,R);

if x=-infinity then if type(y,realcons) then return cat(t, "<=", y) else

cat(t, "<", op(1,y)) fi else

if type(x,realcons) then return cat(t, ">=", x) else

cat(t, ">", op(1,x)) fi; fi;

end proc:

Examples of use:

ConvIntoIneq(not x<100), ConvIntoIneq(not x>100), ConvIntoIneq(not x>=100), ConvIntoIneq(not x<=100), ConvIntoIneq(x>100), ConvIntoIneq(x<100), ConvIntoIneq(x>=100) ;

 Addition. Unfortunately, if we want to get a "nice" output, the name of variable can appear on the right side of the inequality:

map(parse, [ConvIntoIneq(not x<100), ConvIntoIneq(not x>100), ConvIntoIneq(not x>=100), ConvIntoIneq(not x<=100), ConvIntoIneq(x>100), ConvIntoIneq(x<100), ConvIntoIneq(x>=100)])[] ;

I added 2 options into your code: axes=normal and orientation=[15,50]  and used user lighting :

 Addition. The plot was saved in Paint as png - file.

Another variant. Additionally I added  color=khaki ,  style=surface , numpoints=20000  and changed  user lighting  :

plot([[0,0.36],[5,0.38],[10,0.3],[15,0.18],[20,0.96],[25,0.18],[30,0.16],[35,0.96]], style=line, thickness=5, color=red, view=[0..35,0..1], axis[1] = [gridlines = [7, thickness = 1, subticks = false, color = grey]], axis[2] = [gridlines = [6, thickness = 1, subticks = false, color = grey]], tickmarks=[[0=1,5=5,10=10,15=50,20=250,25=1000,30=10000,35=50000], spacing(0.2,0)], size=[650,350], font=[COURIER,ROMAN,14]);

 Edited.

P := (L, a, b) ->plot(L, x = a .. b, color = [seq(`if`(i::odd, red, black), i = 1 .. nops(L))]):   # L is a list of functions

Example of use:

P([sin(x), cos(x), 2, 3-x], 0, 5);

First, you can find the indices of columns that must be deleted, and then delete them.

Example:

A:=<0,1,0,3; 0,-2,0,4>;

ListTools[SearchAll](0, convert(A[1], list));

LinearAlgebra[DeleteColumn](A, [%]);

Your syntax is wrong there. Below is the text of the procedure called  P  that solves the problem for any  Time  and  solution:

P:=proc(Time, solution)

local n;

n:=nops(Time);

piecewise( seq( op([t>=Time[i] and t<Time[i+1], solution[i]]), i = 1..n-1));

end proc:

Solution of above example:

P([0,1,2,3], [t^2, 1, 3-t]);

plot(%, t=0..3, scaling=constrained);

Edited.

Example:

F:=piecewise(t>=0 and t<1, t^2, t>=1 and t<2, 1, t>=2 and t<3, 3-t);

plot(F, t=0..3, scaling=constrained);

Addition:  outside the interval t=0..3  the function  F by default considered to be equal to 0 .