How to write a function in the following form f[a_, b_, c_,...] so that it is flexible with respect to the number of arguments?

Okay, this will be long and involved. I'm basing this as much as I can on your code so that you don't have to try to translate too much. But I've only gotten halfway to your final objective. I hope you can follow this pattern to get the rest of the way (ask further questions if you can't). There are many other improvements we could make, but I wanted to stay focused as much as possible on how to write an extensible function.

Our first objective is to rewrite this:

Kxx[Sequence @@ var1] = Kx 

I'm deferring tackling the Pxx and EER and so forth.

After inspecting your var1, we see that has four arguments at the front and then a variable number of arguments after that. But those variable arguments can always be made with a Table and iteration. That means they're really just a matrix. So, let's try that:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] := Kx 

Note the :=. That suspends evaluation until we need it. Now let's plug in your Kx:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] := -1/2 2 Pi Sum[ c[[i1, j1]] c[[i2, j2]] Kk1, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}] 

Notice I changed [] to [[]], because now we're accessing a matrix. Now consider that imax and jmax can be derived from c. Let's use that fact:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] := With[ {dims = Dimensions[c]}, -1/2 2 Pi Sum[ c[[i1, j1]] c[[i2, j2]] Kk1, {i1, 1, dims[[1]]}, {i2, 1, dims[[1]]}, {j1, 1, dims[[2]]}, {j2, 1, dims[[2]]}]] 

Time to plug in your Kk1, but first let's inspect it on its own:

Integrate[Kk r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}] 

We need definitions for a and b. Your definitions look like this:

a[1] = a1; b[1] = b1; Do[a[i] = a[i - 1] qa, {i, 2, imax}]; Do[b[j] = b[j - 1] qb, {j, 2, jmax}]; 

There is a problem here and an opportunity. The problem is that these definitions refer to a1, qa, etc. But as you see above we're passing in those arguments to Kxx, so we need to feed them to a and b. The opportunity is that a and b are simple functions, so we don't need to use a loop to pre-define values for them. In fact, they're the same function. As you noted, they are geometric sequences. So, let's do this:

geoseq[init_, r_, n_Integer?Positive] := init*r^(n - 1) 

Let's plug this into our Kxx, but let's keep it slightly simpler by redefining a and b locally:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] := With[ {dims = Dimensions[c], a = geoseq[a1, qa, #] &, b = geoseq[b1, qb, #] &}, -1/2 2 Pi Sum[ c[[i1, j1]] c[[i2, j2]] Integrate[Kk r, {r, 0, Infinity}, {z, -Infinity, Infinity}, Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}], {i1, 1, dims[[1]]}, {i2, 1, dims[[1]]}, {j1, 1, dims[[2]]}, {j2, 1, dims[[2]]}]] 

Okay, now for Kk. It looks like this:

FullSimplify[Psi[r, z, i2, j2]*Laplacian[Psi[r, z, i1, j1], {r, \[Theta], z}, "Cylindrical"]] 

So, we need Psi. Your Psi depended on a and b. But we've replaced those with geoseq. And now this should be déjà vu... we need to add arguments:

Psi[a1_, b1_, qa_, qb_, r_, z_, i_, j_] := Exp[-geoseq[b1, qb, j]*z^2]*Exp[-geoseq[a1, qa, i]*r^2] 

Okay, putting it all together, and hoping I can do that without a copy-paste error:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] := With[ {dims = Dimensions[c], a = geoseq[a1, qa, #] &, b = geoseq[b1, qb, #] &}, -1/2 2 Pi Sum[ c[[i1, j1]] c[[i2, j2]] Integrate[ FullSimplify[ Psi[a1, b1, qa, qb, r, z, i2, j2]* Laplacian[ Psi[a1, b1, qa, qb, r, z, i1, j1], {r, \[Theta], z}, "Cylindrical"]] r, {r, 0, Infinity}, {z, -Infinity, Infinity}, Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}], {i1, 1, dims[[1]]}, {i2, 1, dims[[1]]}, {j1, 1, dims[[2]]}, {j2, 1, dims[[2]]}]] 

Let's take it for a spin:

Kxx[1.22, 0.44, 1.40, 1.81, {{1, 2}}] (* 27.5497 *) 

Notice the last argument is a matrix, we needed to reformulate from what you had. Let's try another:

Kxx[1.22, 0.44, 1.40, 1.81, {{1, 2, 1}, {3, 2, 1}}] (* 281.248 *) 

At this point there are lots of opportunities for clean up and maybe some performance improvement, but I'm going to stop there for now.