Basically it's clear that diff doesn't understand dummy variables in sums and products.  On the other hand, if what you have involves finite sums and/or products with explicit bounds, e,g, product(x[i], i=1..3), Maple will expand those out, in this case to x[1]*x[2]*x[3].  The derivatives of an explicit expression such as this should be fine.
My suggestion is to stick with such expressions (or at least test out your results by trying them on such expressions).

Basically it's clear that diff doesn't understand dummy variables in sums and products.  On the other hand, if what you have involves finite sums and/or products with explicit bounds, e,g, product(x[i], i=1..3), Maple will expand those out, in this case to x[1]*x[2]*x[3].  The derivatives of an explicit expression such as this should be fine.
My suggestion is to stick with such expressions (or at least test out your results by trying them on such expressions).

Are you sure your solution "by hand" is correct?  I very much doubt that this could be put into the form exp(-cx)/x.

Are you sure your solution "by hand" is correct?  I very much doubt that this could be put into the form exp(-cx)/x.

Actually there's a more serious bug here.  Since i is a dummy variable in the product, it doesn't make sense to differentiate with respect to z[i].  But the answer Maple returns makes even less sense.  It would be the sum of the derivatives with respect to all the z[i].  Similarly:

> S := sum(z[i], i = 1 .. n);
   diff(S, z[i]);

n

On the other hand:

> diff(S,z[1]);

0

Actually there's a more serious bug here.  Since i is a dummy variable in the product, it doesn't make sense to differentiate with respect to z[i].  But the answer Maple returns makes even less sense.  It would be the sum of the derivatives with respect to all the z[i].  Similarly:

> S := sum(z[i], i = 1 .. n);
   diff(S, z[i]);

n

On the other hand:

> diff(S,z[1]);

0

It's almost certainly impossible to solve an equation such as this in closed form.  You could try numerical solutions (for given values of the parameters), or perhaps series solutions in powers of some small parameter.

It's almost certainly impossible to solve an equation such as this in closed form.  You could try numerical solutions (for given values of the parameters), or perhaps series solutions in powers of some small parameter.

How did it not work?

Note that I omitted the initial definition of A, which was

> A := Matrix(5,5,f);

as in your original code.

This latest way should work if instead of Matrix(5,5,f) you used

Matrix(5,5,(i,j)->f(i,j));

The reason this works but Matrix(5,5,f) does not is rather subtle, I think.

Initially f was undefined, but after assigning values to its remember table it does have a value:

> eval(f);

   proc() option remember; 'procname(args)'  end proc;

 My guess is that when you call Matrix(5, 5, f), what actually gets applied to make each matrix element is not f itself but its value, no longer attached to the name f.
And then the procname is unknown rather than f. 

How did it not work?

Note that I omitted the initial definition of A, which was

> A := Matrix(5,5,f);

as in your original code.

This latest way should work if instead of Matrix(5,5,f) you used

Matrix(5,5,(i,j)->f(i,j));

The reason this works but Matrix(5,5,f) does not is rather subtle, I think.

Initially f was undefined, but after assigning values to its remember table it does have a value:

> eval(f);

   proc() option remember; 'procname(args)'  end proc;

 My guess is that when you call Matrix(5, 5, f), what actually gets applied to make each matrix element is not f itself but its value, no longer attached to the name f.
And then the procname is unknown rather than f. 

You might try something like this.

> X := proc ... end proc:
  XInert:= ToInert(eval(X));
  params:= {op(op(1,XInert))};
  locals:= {op(op(2,XInert))};
  globals:= indets(XInert,_Inert_NAME(anything)) 
       minus params minus locals; 
  newlocals:= locals union 
     remove(t -> type(convert(op(t),name),protected), globals);
  NewX:= FromInert(subsop(2=_Inert_LOCALSEQ(op(newlocals)),XInert));
  

Note: by "all" I assume you don't really mean "all names used in the procedure are local".  For example, if your procedure uses the sin function you don't want sin to be a local variable.  So in the above I just made every non-protected name local.

You might try something like this.

> X := proc ... end proc:
  XInert:= ToInert(eval(X));
  params:= {op(op(1,XInert))};
  locals:= {op(op(2,XInert))};
  globals:= indets(XInert,_Inert_NAME(anything)) 
       minus params minus locals; 
  newlocals:= locals union 
     remove(t -> type(convert(op(t),name),protected), globals);
  NewX:= FromInert(subsop(2=_Inert_LOCALSEQ(op(newlocals)),XInert));
  

Note: by "all" I assume you don't really mean "all names used in the procedure are local".  For example, if your procedure uses the sin function you don't want sin to be a local variable.  So in the above I just made every non-protected name local.

OK, you want an antiderivative.  No, you don't need to define limits.

>  int(exp(-c*x)/sqrt(a^2+b*x^2), x);

int(exp(-c*x)/sqrt(a^2+b*x^2), x)

Maple is telling you it can't find a formula for the antiderivative.  I also tried it on integrals.com, and that didn't find one either.  Most likely there is no closed-form formula.

OK, you want an antiderivative.  No, you don't need to define limits.

>  int(exp(-c*x)/sqrt(a^2+b*x^2), x);

int(exp(-c*x)/sqrt(a^2+b*x^2), x)

Maple is telling you it can't find a formula for the antiderivative.  I also tried it on integrals.com, and that didn't find one either.  Most likely there is no closed-form formula.

This code is iterative, not recursive.

A recursive function would call itself (with different inputs)

An iterative function contains a loop: in this case

while Err > s do ... end do