In the procedure  eul   in these lines

  x:= 0 + j*h:

  y := y[j]:

you assign new values to the formal parameters  x  and  y . It is invalid.

In your equation all the constants should be exact (integers, fractions and so on) rather than floats.

Two examples:

solve(3/5*x+0.1=0);

solve(3/5*x+1/10=0);

This is a known problem of mathematical folklore! I saw it in different options (beetles, flies, turtles, etc.), sometimes a regular triangle occurs instead of a square.

The basic ideas of the solution are: 1) the use of symmetry, on the basis of which the meeting is always in the center of the corresponding regular polygon and 2) velocity decomposition in 2 components, one goes to the center, and the second is perpendicular the first.

Pretty simple task there finding the time before the meeting - you just need to divide the distance to the center by the corresponding component of the velocity. In this example, we get  t = (L/sqrt(2))/(v/sqrt(2))=L/v  (L is the edge of the square and  v is the velocity of a dancer). Now it is easy to find the length of any trajectory by multiplying the speed of the dancer by the time. We get  S=L (in the example  L=4). 
Much more interesting, of course, set the trajectory of the dancer, and to find the laws of motion, that is, the function determines the position of the dancer on plane at any time. To find the equations of the trajectories necessary to solve differential equations (fairly simple equations with separable variables). The trajectories of the dancers themselves do not depend on the speed v,  only the time of movement along these trajectories does. These paths are called logarithmic spirals.

Here is solution with Maple together with animation:

restart;

for i to 4 do

sol[i] := dsolve({r(-2*Pi*(i-1)*(1/4)) = L/sqrt(2), D(r)(phi) = r(phi)}, r(phi));

od:

seq(sol[i], i = 1 .. 4);  # Equations of the trajectories in polar coordinates (phi<=0) 

assign(seq(r[i] = unapply(rhs(sol[i]), phi), i = 1 .. 4));

simplify(int(sqrt(r[1](phi)^2+(D(r[1]))(phi)^2), phi = -infinity .. 0), assume = positive); # Finding of the length of one trajectory

for i to 4 do

Sol[i] := dsolve({r[i](phi(t))*(D(phi))(t) = -v/sqrt(2), phi(0) = -2*Pi*(i-1)*(1/4)}, phi(t)):

od:

seq(Sol[i], i = 1 .. 4); # The laws of changing  phi(t)  ( t  is the time)

assign(seq(phi[i] = unapply(rhs(Sol[i]), t), i = 1 .. 4));

# The animation 

L := 4: v := 1: L1 := 4.2:

for i to 4 do

Room := plottools[curve]([[L1/sqrt(2), 0], [0, L1/sqrt(2)], [-L1/sqrt(2), 0], [0, -L1/sqrt(2)], [L1/sqrt(2), 0]], color = blue, thickness = 3), plottools[polygon]([[L1/sqrt(2), 0], [0, L1/sqrt(2)], [-L1/sqrt(2), 0], [0, -L1/sqrt(2)]], color = "LightYellow");

A[i] := plots[animate](plot, [[r[i](phi[i](t))*cos(phi[i](t)), r[i](phi[i](t))*sin(phi[i](t)), t = 0 .. T], color = green, thickness = 3], T = 0 .. L/v, scaling = constrained, frames = 80);

B := (x, y)->plots[pointplot]([[x, y]], color = red, symbol = solidcircle, symbolsize = 25);

C[i] := plots[animate](B, [r[i](phi[i](T))*cos(phi[i](T)), r[i](phi[i](T))*sin(phi[i](T))], T = 0 .. L/v, frames = 80)

od:

plots[display](Room, A[1], A[2], A[3], A[4], C[1], C[2], C[3], C[4]);

 Dancers.mw

f:=x^2*exp(-1.2*x);

maximize(f, x=0..100, location);

a:=op([2,1,1,1,2], [%]);

For exact solution replace  1.2  with  6/5

Addition:  another variant

f:=x^2*exp(-6/5*x);

maximize(f, x=0..100, location);

a:=rhs(%[2,1,1,1]);

If you have any parameters it is necessary to specify the unknown function. 

Examples:

dsolve({diff(f(x),x,x)-3*diff(f(x),x)+2=x^2+1, D(f)(a)=b}, f(x));

dsolve({diff(f(x),x,x)-3*diff(f(x),x)+2=x^2+1, D(f)(a)=b});

Since your equation is of second order must be another initial condition.

Example of solution with the plot:

R := 5;

g := 9.8;

deq1 := {diff(x(t),t,t) = -g*sin(x(t))/R, x(0) = 20, D(x)(0)=2};

dsol1 := dsolve(deq1, numeric);

plots[odeplot](dsol1,[t,x(t)], t=0..20, color=red, thickness=2);

I think you are interested in non trivial (non-zero) solutions. So use  solve  instead of  fsolve . solve command finds all the solutions of a polynomial system in the form  RootOf. Using then  evalf  command, you can numerically to find these solutions.

Example:

restart;

A, U, U1:=LinearAlgebra[RandomMatrix](4,4,generator=-9..9), <seq(u[i], i=1..4)>, <seq(u[i]^2, i=1..4)>;

sys:=Equate(A.U1, U);

sol:=solve(%):

fsol:=evalf(%);

eval((lhs-rhs)~(sys), fsol[2]);  # Verification

Simple procedure  cplot  generalizes  Acer's way and shows also by specific color the points and ranges, in which the function is equal to zero. The colors and thickness of the points, lines and sizes of the points are set by default, but can also be specified by the user.

cplot := proc (f::algebraic, r::range, Colors::list := [blue, red, green], R::realcons := 17, t::`=` := thickness = 2)

local x, L, M, M1, Points, P;

uses plots;

x := indets(f)[1];

L := map(op, {solve({lhs(r) <= x, x <= rhs(r), f}, allsolutions, explicit)});

M := select(s->type(s, `=`) or type(s, `<=`), L);

M1 := select(t->type(t, realcons), map(op, M));

Points := plot([seq([M1[i], 0], i = 1 .. nops(M1))], style = point, symbol = solidcircle, symbolsize = R, color = Colors[3]);

display(plot([piecewise(0 < f, f, undefined), piecewise(f < 0, f, undefined), piecewise(f = 0, f, undefined)], x = r, color = Colors, t), Points);

end proc:

Example:

F := (x+1)*x*(x-2):

f := piecewise(x > -1 and x < 0, 0, F);

cplot(f, -2 .. 3, thickness = 3);

The procedure  cplot  is easily modified to show the intervals of increase, decrease, and constancy of the function.

cplot1 := proc (f::algebraic, r::range, Colors::list := [blue, red, green], R::realcons := 17, t::`=` := thickness = 2)

local x, f1, L, M, M1, Points, P;

uses plots;

x := indets(f)[1];

f1 := diff(f, x);

L := map(op, {solve({lhs(r) <= x, x <= rhs(r), f1}, allsolutions, explicit)});

M := select(s->type(s, `=`) or type(s, `<=`), L);

M1 := select(t->type(t, realcons), map(op, M));

Points := plot([seq([M1[i], eval(f, x = M1[i])], i = 1 .. nops(M1))], style = point, symbol = solidcircle, symbolsize = R, color = Colors[3]);

display(plot([piecewise(0 < f1, f, undefined), piecewise(f1 < 0, f, undefined), piecewise(f1 = 0, f, undefined)], x = r, color = Colors, t), Points);  

end proc:

Example:

F := (x+1)*x*(x-2) - 1:

f := piecewise(x > -1 and x < 0, -1, F);

cplot1(f, -2 .. 3);    # Green color - the points in which the derivative is equal to zero

Another example:

f := 2*sin(t)+3*cos(t)-1:

cplot1(f, -Pi .. 2*Pi, thickness = 3) ;

Edit. Codes have been edited and added another example.

For presentation the different surfaces are painted in different colors:

restart;

A := plots[spacecurve]([x, 0, ln(x)], x = 1 .. exp(1), thickness = 2, color = red):

B := plots[spacecurve]([x, 0, -ln(x)], x = 1 .. exp(1), thickness = 2, color = blue):

C := plots[spacecurve]([exp(1), 0, z], z = -1 .. 1, thickness = 2, color = yellow):

N := 100:

S1 := seq(plot3d([[x*cos(phi), x*sin(phi), ln(x)], [x*cos(phi), x*sin(phi), -ln(x)]], x = 1 .. exp(1), phi = 0 .. 2*Pi*k/N, color = [red, blue], style=surface), k = 1 .. N):

S2 := seq(plot3d([exp(1)*cos(phi), exp(1)*sin(phi), z], phi = 0 .. 2*Pi*k/N, z = -1 .. 1, color = yellow, style=surface), k = 1 .. N):

S := plots[display](A, B, C), seq(plots[display](A, B, C, S1[k], S2[k]), k = 1 .. N):

plots[display](S, insequence=true, scaling = constrained, view = [-4.1 .. 4.1, -4.1 .. 4.1, -1.4 .. 1.4], orientation = [-75, 60], axes = normal);

Your surface is composed of two parts: external and internal. The outer part is indeed a sphere and it covers the inner surface. I made a cut-out and painted in different colors, so that you can see the inner surface. Your body is located between these surfaces.

A := plot3d([x, (-x^2+2*x)*cos(phi), (-x^2+2*x)*sin(phi)], x = 0 .. 2, phi = -Pi .. (1/2)*Pi, color = yellow):

B := plot3d([2*cos(alpha)^2, 2*cos(alpha)*sin(alpha)*cos(phi), 2*cos(alpha)*sin(alpha)*sin(phi)], alpha = 0 .. (1/2)*Pi, phi = -Pi .. (1/2)*Pi, color = cyan):

plots[display](A, B, axes = normal, lightmodel = light4, orientation = [-115, 55]);

Procedures  Y_low, Y_up, Y for each value of  X  are respectively looking  for  Y  for the lower branch, the upper branch and both values. The procedures are based on successive elimination of the variable  theta, beginning with the first equation without  elimination  command. I tried to write everything in one procedure, but problems arose when plotting.

restart;

f := theta->cos(theta)+(2/25)*cos(3*theta):

g := theta->-sin(theta)+(2/25)*sin(3*theta):

Y_low := proc (X)

local m, M, c;

m := minimize(f(theta)); M := maximize(f(theta));

if is(X < m) or is(M < X) then return "No solutions" end if;

c := RealDomain[solve](subs(cos(theta) = y, expand(f(theta) = X)));

simplify(subs({cos(theta) = c, sin(theta) = sqrt(1-c^2)}, expand(g(theta))));

end proc:

Y_up := proc (X)

local m, M, c;

m := minimize(f(theta)); M := maximize(f(theta));

if is(X < m) or is(M < X) then return "No solutions" end if;

c := RealDomain[solve](subs(cos(theta) = y, expand(f(theta) = X)));

simplify(subs({cos(theta) = c, sin(theta) = -sqrt(1-c^2)}, expand(g(theta))));

end proc:

Y:=X->[Y_low(X), Y_up(X)]: 

Examples:

Y_low(1); evalf(%);

Y_up(1); evalf(%);

Y(0);

Y(2);

plots[display](plot(Y_up, -0.5 .. 1, color=red, thickness=2, scaling=constrained), plot(Y_low, -0.5 .. 27/25, color=blue, thickness=2, scaling=constrained));

with(plots):

A := spacecurve([-y^2+4, y, 2-y], y = 0 .. 2, color = red, thickness = 3):

B := plot3d(2-y, x = 0 .. 4, y = 0 .. 2, color = green, style=surface):

C := plot3d([-y^2+4, y, z], y = 0 .. 2, z = 0 .. 2, color = cyan, style=surface):

display(A, B, C, scaling = constrained, axes = normal, view = [-1 .. 4.5, -1 .. 2.4, -1 .. 2.4], lightmodel=light4);

 Addition.  The easiest way of cutting away - to use  filled  option:

restart;

A := plots[spacecurve]([-y^2+4, y, 2-y], y = 0 .. 2, color = red, thickness = 3):

B:=plot3d(2-y, x = 0..4-y^2, y=0..2, color = grey, style=surface, filled, scaling=constrained, lightmodel=light4, axes=normal, view = [-1 .. 4.5, -1 .. 2.4, -1 .. 2.4]):

plots[display](A,B); 

restart;

f :=(theta)->cos(theta)+0.8e-1*cos(3*theta);

g :=(theta)->-sin(theta)+0.8e-1*sin(3*theta);

X^2=f(theta)^2;

Y^2=g(theta)^2;

simplify(expand(%%+%));

solve(%, Y);

You have one nonlinear transcendental equation with 6 unknowns. In general, the set of solutions of similar equations - nonlinear 5-dimensional manifold in the 6-dimensional space. It may turn out to be empty. But in your example, there are infinitely many solutions. The variable  a2  will be regarded as the main variable, other variables will be considered as parameters. Since a2 is included in the equation by linear way, then for each permissible  set of parameters exists not greater than one solution of the equation with respect to a2 matched your criteria. Here is an example of finding one solution (for exact solution all constants should be exact):

restart;

Eq:=50*tan(8*Pi/180)=(a1/(d2-d1))*ln(d2/d1)+(a2/(d3-d2))*ln(d3/d2)+((6/10-a1-a2)/(d4-d3))*ln(d4/d3);

solve(eval(Eq,[d1=5/100, d2=6/100, d3=7/100, d4=9/10, a1=1/10]));

evalf(%);

If you need more solutions, the main problem is the selection of parameters for which and around which there are valid solutions. This can be done in the usual for loop:

N:=0:

for i from 0.01 to 0.06 by 0.01 do

for j from i+0.01 to 1 by 0.1 do

for k from j+0.01 to 1.5 by 0.1 do

for m from k+0.01 to 2 by 0.1 do

for n from 0.1 to 0.5 do

p:=fsolve(eval(Eq,{d1=i,d2=j,d3=k,d4=m,a1=n}));

if p>0 and n+p<0.6 then N:=N+1; L[N]:=[i,j,k,m,n,p] fi;

od: od: od: od: od:

L:=convert(L,list);

Easy modification of  coefff  procedure solves the problem:

coefff:=proc(P, T, t)

local L, H, i, k, s, t1:

L:=[coeffs(P, T, 'h')]: H:=[h]: k:=0:

s:=`if`(t::realcons,t,op(1,t));

if s::realcons and s<>1 then t1:=t/s else t1:=t fi;

for i from 1 to nops(H) do

if H[i]=t1 then k:=L[i] fi:

od:

if t1=t then k else k/s fi;

end proc:

Examples:

f:=a*x^2+b*x*y^3+2;

coefff(f,[x,y],x^2/2);

coefff(f,[x,y],x*y^3);

coefff(f,[x,y],1/10);