@Carl Love Thank you for the response! Your method is a much better way than using string tools. I can get this to work for simple functions like z*x shown in my original post. However, the actual function I want to find the matrix output of values for is much longer with many elliptical integrals. For some reason this method will not work to evaluate it at different values of z and x. Instead of a numerical output, I get a matrix that still contains the variables z and x. Maybe the output is not evaluated due to some conditional expressions where z and x must be within certain bounds.

Do you know why this could be? Attached is my formula for what I actually want to calculate (also changing x to r). 

isol:=10;
Zres:=3:
Rres:=2:
l:=1;
d:=0.25;
mu[0]:=4*Pi*(10^(-7)):
R:=d/2:
ires:=(isol*Nturns*Nsets)/Zres:
kk:=(4*R*r)/(((R+r)^2)+((z-A)^2)):
K:=int(1/(sqrt((1-(t^2))*(1-(kk)*(t^2)))),t=0..1):
E:=int(sqrt(1-(kk)*(t^2))/(sqrt(1-(t^2))),t=0..1):
Bz:=((mu[0]*ires)/(2*Pi))*(1/(sqrt(((R+r)^2)+((z-A)^2))))*(K+E*((R^2)-(r^2)-((z-A)^2))/(((R-r)^2)+((z-A)^2))):
BzL:=eval~(Bz,A=~[seq((n*l)/(Zres),n=((-Zres)/2)..((Zres)/2))]):
Bzsum:=(sum(BzL[n],n=1..(Zres))):
with(LinearAlgebra):
M:=(z,r)->Bzsum;
M1:=Array(1..Zres,i->(i*l)/(2*Zres));
M2:= Array(1..Rres,i->(i*d)/Rres);
A=Array(1..numelems(M1),1..numelems(M1),(i,j)->M(M1[i],M2[j]));