There are several things going on with your existing code. The first tip is that using Subscript can cause complications, so people tend to recommend that beginners don't use it for variable names.

The second tip is that the VectorAnalysis package is obsolete since version 9, and in any case, the usual Dot (.) function built into the main system will do what you want.

The third (and more fundamental) issue is that you have not actually defined functions of x and y for your DensityPlot to address. The x and y used in a plotting function are not the same as the ones you use outside that plot to define the expression to be plotted. This is a scoping issue; I recommend you have a look at some of the documentation on this issue. You instead need to define a function that takes a parameter, using the SetDelayed notation (:=). You would, for example, rewrite your first $\psi$ function as:

psis[x_] := 
 Cos[DotProduct[k1, r]]*x + Cos[DotProduct[k2, r]]*x + Cos[DotProduct[k3, r]]*x;

Or using the built-in Dot:

psis1[x_] := Cos[k1.r]*x + Cos[k2.r]*x + Cos[k3.r]*x

Re-writing the function shows another issue: your functions are linear in x and y respectively. I think you actually want to rewrite them as follows.

psis2[x_] := Cos[k1.r*x] + Cos[k2.r*x] + Cos[k3.r*x]; 
psia2[y_] := Sin[k1.r*y] + Sin[k2.r*y] + Sin[k3.r*y];

This becomes obvious when you plot each version (note I've added something to psis1 because it's otherwise identical to psis):

Plot[{psis[x], psis1[x] + 10, 5 psis2[x]}, {x, 0, 100}]

So using these changed versions of your functions, we get closer to what you expect:

DensityPlot[(psis2[x])^2 + ( psia2[y])^2, {x, 0, 100}, {y, 0, 100}, 
  ColorFunction -> GrayLevel]

It's not quite as in your desired picture, which leads me to suspect that you need to tweak your functional forms a bit more. But as I'm an economists and not a physicist, I'll have to leave that bit to you.

There are several things going on with your existing code. The first tip is that using Subscript can cause complications, so people tend to recommend that beginners don't use it for variable names.

The second tip is that the VectorAnalysis package is obsolete since version 9, and in any case, the usual Dot (.) function built into the main system will do what you want.

The third (and more fundamental) issue is that you have not actually defined functions of x and y for your DensityPlot to address. The x and y used in a plotting function are not the same as the ones you use outside that plot to define the expression to be plotted. This is a scoping issue; I recommend you have a look at some of the documentation on this issue. You instead need to define a function that takes a parameter, using the SetDelayed notation (:=). You would, for example, rewrite your first $\psi$ function as:

psis[x_] := 
 Cos[DotProduct[k1, r]]*x + Cos[DotProduct[k2, r]]*x + Cos[DotProduct[k3, r]]*x;

Or using the built-in Dot:

psis1[x_] := Cos[k1.r]*x + Cos[k2.r]*x + Cos[k3.r]*x

Re-writing the function shows another issue: your functions are linear in x and y respectively. I think you actually want to rewrite them as follows.

psis2[x_] := Cos[k1.r*x] + Cos[k2.r*x] + Cos[k3.r*x]; 
psia2[y_] := Sin[k1.r*y] + Sin[k2.r*y] + Sin[k3.r*y];

This becomes obvious when you plot each version (note I've added something to psis1 because it's otherwise identical to psis):

Plot[{psis[x], psis1[x] + 10, 5 psis2[x]}, {x, 0, 100}]

So using these changed versions of your functions, we get closer to what you expect:

DensityPlot[(psis2[x])^2 + ( psia2[y])^2, {x, 0, 100}, {y, 0, 100}, 
  ColorFunction -> GrayLevel]

It's not quite as in your desired picture, which leads me to suspect that you need to tweak your functional forms a bit more. But as I'm an economist and not a physicist, I'll have to leave that bit to you.