Numerical solution for a system of linear equations

Hello,

I have a system of linear equations which I am trying to solve numerically. The physical problem I am trying to solve is a Laplace equation in a specific topology (a rectangle with some perks in it) but I think the problem is more general. The solution is the list of coefficients for the series composing the solution for the Laplace equation.

The way I solve: The equations I am trying to solve are of the form: sum(A_n*Sinh(x_i a_n)*sin(z1*a_n), n)=b, when I am solving for A_n. The x_i are points I choose from a range of [0,r]. I can choose an arbitrary number of points, so I can have an arbitrary number of equations. When I choose n equations for n variables, I get a solution, but it changes dramatically with the choice of x_i (it really shouldn't! it should give the same values no matter the choice of x_i). When I add more equations, for the same n variables, I don't get a solution.

It is really crucial for me to solve these equations so every help will be much appreciated.

Here is the code:

%dimensions of the box
r=0.4;
R=1;
z1=0.6;
z2=1;
zMax=1;
%the voltages on the edges
V0=100;
V1=100;
V2=100;
V3=100;
V4=100;
V5=100;
n=70;%number of variables to solve for
nVec=1:n;
eqNo=n*2;%number of equations

%defining the symbolic variables 
vars=[sym('C1_',size(nVec)) sym('A3_',size(nVec))];
z=sym('z');
bnVec=nVec*pi/z1;
anVec=nVec*pi/r;
cnVec=nVec*pi/(R-r);
zVals=linspace(0.01,z1*0.99,eqNo-n);
%a vector which sets to zero the contribution of the even coefficients
zeroEven=ones(size(nVec));
zeroEven(2:2:n)=0;
%define the equations
eqns=sym(zeros(eqNo,1));
%create the potential equations
for i=1:n
    eqns(i)=vars(i)*sinh(bnVec(i)*r)-vars(i+n)*sinh(bnVec(i)*(R-r))==0;
end
%create the derivative equations
coeffVec(z)=sum(bnVec.*cosh(bnVec*r).*sin(bnVec*z).*vars(1:n)+bnVec.*cosh(bnVec*(R-r)).*sin(bnVec*z).*vars(n+1:2*n));%+A1*(z+C1)-A3*(z+C3);
 
b(z)=sum((4*V5./(z1*bnVec.*sinh(bnVec*r)).*bnVec.*sin(bnVec*z)+4*V0./(r*anVec.*sinh(anVec*z1)).*anVec.*sinh(anVec.*(z1-z))+4*V1./(r*anVec.*sinh(anVec*z1)).*anVec.*sinh(anVec*z)+...
    4*V0./((R-r)*cnVec.*sinh(cnVec*z1)).*cnVec.*sinh(cnVec*(z1-z))+4*V3./(z1*bnVec.*sinh(bnVec*r)).*bnVec.*sin(bnVec*z)+4*V4./((R-r)*cnVec.*sinh(cnVec*z1)).*cnVec.*sinh(cnVec*z)).*zeroEven);

for i=1:size(zVals,2)
    eqns(n+i)=coeffVec(zVals(i))-b(zVals(i))==0;
end
a=vpasolve(eqns,vars);
%tried linsolve but it didn't work any better
% [A,B] = equationsToMatrix(eqns,vars);
% X = linsolve(A,B);
% X=A\B;