Here's your program---debugged, cleaned up, and modernized:

Download NevilleInterp.mw

otat:= (a::`=`, b::list(name=algebraic), c)-> map2(eval, isolate(a,y), b):

otat(x-y = 1, [x=10, x=11]);

     [y = 9, y = 10]

I've translated the algorithm that you posted into Maple. However, the algorithm does not seem to be correct. Here's the translation:

`mod/SFF`:= proc(f::polynom, p::posint)
local
     i:= 1,
     R:= 1,
     x:= indets(f, name)[],
     c, w, g, y, z
;
     if x=NULL then return f mod p end if;
     g:= diff(f,x);
     c:= Gcd(f,g) mod p;
     w:= Quo(f,c,x) mod p;
     while w <> 1 do
          y:= Gcd(w,c) mod p;
          z:= Quo(w,y,x) mod p;
          R:= R*z^i;
          i:= i+1;
          w:= y;
          c:= Quo(c,y,x) mod p
     end do;
     R * (c mod p)
end proc:

The usage is like this:

f:= x^6 + randpoly(x):

SFF(f) mod 3;

To compare with Maple's prepackaged algorithm, use

Sqrfree(f) mod 3;

Your worksheet prompts the user for input. That is a very, very crude way (IMO) to use/write a Maple program. IMO, the only acceptable way for a Maple program to accept input is through a procedure's parameters.

Your wa and wb should be w*a and w*b.

When setting the values of the parameters a, b, and w, you need to use := rather than =.

On the right side of eqx, you use "naked" x; this should be x(t).

On the right side of eqz, you use y; I'm guessing that this should be z(t).

The squaring operation can also be done on a list of lists using elementwise operators:

(`^`~)~(L,2);

For editing code in a worksheet, try Code Edit Regions. You can find them on the Insert menu. 

I see that you're using Maple 18. There are two nuisances with Code Edit Regions in Maple 18. I think that the first of these has been corrected in Maple 2015 and that the second hasn't. The first nuisance is that when you execute the code in the Region, if there's an error, it doesn't tell you where the error is nor does it place the cursor near the error. The second nuisance is that like most worksheet-embedded components, it's very hard to get your mouse cursor precisely on the control points that let you resize the Region, which you'll very likely need to do.

I don't know about Volterra equations, but your integro-differential equation can be easily converted to a third-order ODE-IVP which can be easily solved symbolically. So, I say, "Why bother with Volterra equations?" To do the solving entirely in Maple:

Download IntegroDiffEqn.mw

Here's another way to do it without unapply, using subs. This method is extremely powerful for three reasons:

1. Whereas unapply can only be applied to a single expression, this subs method can be applied to any procedure whatsoever.
2. This method can be used after the fact: i.e., after f, g, and h have already been defined.
3. It is very fast, based on the built-in command subs.

However, my personal opinion is that imprudent use of this method can make code difficult to read. So, for this particular example, I'd use Kitonum's second method (my first choice) or unapply.

The subs method is shown in the line offset by (**) below.

f:= 2*x:  g:= x^2:
h:= x-> piecewise(x <= 5, _f, _g):
(**) h1:= subs([_f= f, _g= g], eval(h)); (**)
h1~([1,2,8]);

I show two other unrelated things in the above code fragment, one related to map and the other to piecewise. The former is that map can often be replaced by the elementwise operator ~. See ?elementwise. The latter is that inequalities compounded with and are usually unnecessary (and wasteful) as piecewise conditions. The conditions are processed left to right; the fact that the nth condition is being evaluated implies that conditions 1, ..., n-1 are false. I very often see even experienced Maple programmers using these unnecessary compound conditions in piecewise. 

Since both the commands D and diff already provide this functionality, I don't know why you'd want to create the procedure MultiDiff. But, assuming that you have some reason for doing it, ...

MultiDiff:= proc(f, n::posint, x::name) (D@@n)(f)(x) end proc;

The use command has several subleties; I find it very difficult to use (pun intended). Use fully qualified procedure names instead. That is, replace references to Color with ColorTools:-Color.

You've made a common Maple error that has nothing to do with either map or piecewise. I don't know if this error has an official name, so I'll call it an attempted indirect reference to a procedure parameter. (Perhaps that'll become its official name.)

Your h is a procedure (also called a function) whose parameter is x. Your expressions f and g depend on global x. It's not the same x. When you placed f and g within h's body, you were hoping that the indirect reference to x would be to the parameter x. It doesn't work that way. The command unapply is used to resolve these indirect references so that all references to x within an expression refer to the parameter x. Unfortunately, the help page ?unapply doesn't clearly state that.

The command that you want is coeff rather than coeffs:

coeff(r, x, 1);

will return the coefficient of the first power of x. You may then pass this to solve.

Here's the much faster way that I promised---the way that doesn't need to trace back over the path to determine that it'll self intersect on the next step. This generates 1024 walks in an eighth of a second on my machine, the memory usage is insignificant, and the code is much simpler. Indeed, there are no if statements and a only single for loop of two statements.

(Upon re-reading your Question, I see that you're using the word "line" to refer to what the rest of the English-speaking mathematical world (AFAIK) refers to as a random walk. To avoid confusion, I use "walk" in the following. You may read this to yourself as "line". The term path is reserved (by the English-speaking mathematical world) for walks that don't intersect themselves.)

restart:

(M,N):= (2^10, 2^7): #(num walks, max length of walk)

#Tom Leslie's "turn" functions, simplified slightly:
gostraight:= (p,q)-> 2*~p-~q:
turnLeft:=     (p,q)-> p+~`if`(p[1]=q[1], [q[2]-p[2], 0], [0, p[1]-q[1]]):
turnRight:=   (p,q)-> p+~`if`(p[1]=q[1], [p[2]-q[2], 0], [0, q[1]-p[1]]):
dirs:= [gostraight, turnLeft, turnRight]:
R3:= rand(1..3):

OneWalk:= proc(MaxLength::posint)
local
     W:= table([1= [0,0], 2= [1,0]]), #The walk
     V:= table([[0,0]= NULL]), #The points visited, without the terminus
     k
;
     for k from 3 to MaxLength+1 while not assigned(V[W[k-1]]) do
          V[W[k-1]]:= NULL; #Mark previous point as visited.
          W[k]:= dirs[R3()](W[k-1], W[k-2]) #Generate next point.
     end do;
     convert(W, list)
end proc:

Walks:= CodeTools:-Usage(['OneWalk(N)' $ M]):

memory used=17.46MiB, alloc change=32.00MiB, cpu time=125.00ms, real time=125.00ms, gc time=0ns

They can be plotted by simply plot(Walks); however, 1024 is probably too many to appreciate on a plot, so you can do, for example, plot(Walks[100..150]);

I've found all the roots in one period of the function in under a second. Since the function is periodic, this is essentially the same as finding all the real roots.

First, I do a little simplification with the combine command.

restart:
Digits:= 15:
f:= 2*cos(0.5*x)*sin(0.5*x)*cos(3.775*x) +
    2.2075*cos(0.5*x)^2*sin(3.775*x) -
    0.453*sin(0.5*x)^2*sin(3.775*x)
:
fc:= convert(combine(f), rational);

Several things are obvious from this simplified form (several of which were also obvious from the unsimplified form):

1. The function is odd (in the symmetry sense), so the negative of every positive root will also be a root.
2. Zero is a root.
3. The function is periodic with period 80*Pi.
4. The function is well-behaved. It is entire (no discontinuities) and finding roots should be easy.
5. Items 1 and 3 imply that finding the roots on 0..40*Pi will be sufficient.

My first choice for finding all the real roots of such a function is RootFinding:-NextZero.

ffc:= unapply(fc,x):
R:= table(): #To store the roots
r:= 0: #First root
st:= time():
while is(r < 40*Pi) do
     R[r]:= [][]; #Right side can be anything.
     r:= RootFinding:-NextZero(ffc, r, guardDigits= 3);
     if not r::numeric then
          print(r);
          break
     end if
end do:
time()-st;

     0.670

R:= {indices(R, 'nolist')}; #Make table into a set.

(... large set of roots omitted...)

nops(R); #Count the roots.

     191

max((abs@ffc)~(R)); #Check residuals

Since I used three guardDigits, these roots are only guaranteed accurate to 12 digits.

R:= evalf[12](R):

There's no single command to do a double break, but it's a fairly common operation. Here's a standard way of handling it:

for i to M do
     broke:= false;
     for j to 40 while not broke do
          R[i,j];
          for k to j-2 do
               if R[i,j]=K[i,k] then
                    x:= j;
                    print(b[i]=j);
                    broke:= true;
                    break
               end if;
               L[i,j]:= R[i,j]
          end do;
          #The following line isn't needed unless there's code in the j
          #loop after the end of the k loop.
          if broke then break end if;
          # ...possible additional j-loop code.
     end do
end do;

For aesthetic reasons, I usually think of some way to avoid the double break, but I can't state a standard way to avoid it in all cases. If aesthetics are of no concern, then I see no reason to not do it. (But see my "Better ways" Reply below.)