Applying a multi qubit quantum gate to specific qubits
Don't use the swap method, it's very inefficient. And the other person's answer is specific to the CNOT gate, and to be frank, over-complicates things.

Here's a very simple algorithm that solves your problem for every single case, not just the CNOT gate, for any arbitrary bits.

The algorithm:

let sys = matrix representing the current state of the system
let n = number of qubits being simulated
let lgm = logic gate matrix of size 2^n by 2^n
let f = our logic gate transformation function
for i = 0 to (2^n) - 1:
    lgm[column = i] = f(i)
sys = sys × lgg


In classical computers, there's something known as the "decoder". Let's say I have only 3 wires, effectively, 3 bits. But I want to control 8 wires. Can that be done? Yes, because 3 bits has 8 different possibilities: 000, 001, 010, 011, 100, 101, 110, 111. So we can assign each possibility to one of our 8 output wires. This is called "decoding".

If I pass in the number 101 and we have 3 bits, we know 101=5, so then I will be setting output wire 5 to a high voltage and the other 7 output wires will be 0, which we can represent like this: decode(101) = [0, 0, 0, 0, 0, 1, 0, 0].

In this algorithm, I mention the "transformation function" called "f". For classical computers, the transformation function is just taking in an input value and returning the "decoded" version of the output value. So if we have 3 bits, and the output is 4, then we will be returning [0, 0, 0, 0, 1, 0, 0, 0]. We then assign that as the column of our matrix for that value.

Let's think of "decoding" in terms of qubits. How can we decode the qubits |101>?

We know that for our qubit probability vectors, |0> is [1, 0] and |1> is [0, 1]. Decoding qubits can then be done what's called the Kronecker product.

So if we convert each bit to the probability vector equivalent and take the Kronecker product of all of them, we get...

|101&gt; = |1&gt; ⊗ |0&gt; ⊗ |1&gt; = [0, 1] ⊗ [1, 0] ⊗ [0, 1] = [0, 0, 0, 0, 0, 1, 0, 0]


This is how we'd want to decode qubits. This algorithm can be applied to qubits just the same using this.

Let's try this algorithm on a simpler problem.

Let's assume we have a system with only 2 qubits. If we want to apply the Hadamard gate to just 1 qubit, we can generate a logic gate for both qubits that only applies the Hadamard gate to a single qubit. Let's assume the single qubit we want to apply it to is our most significant qubit and the least significant will be unaffected. 

We want a transformation function that for each of our possible inputs, produces the correct output. We have two qubits, this means there are four possible outputs.

f(|00&gt;) = ?
f(|01&gt;) = ?
f(|10&gt;) = ?
f(|11&gt;) = ?


We know the least significant qubit will be unaffected, so we can fill that in.

f(|00&gt;) = ? ⊗ |0&gt;
f(|01&gt;) = ? ⊗ |1&gt;
f(|10&gt;) = ? ⊗ |0&gt;
f(|11&gt;) = ? ⊗ |1&gt;


We also know what the Hadamard does to a qubit, such that:

H(|0&gt;) = 1/sqrt(2)|0&gt; + 1/sqrt(2)|1&gt;
H(|1&gt;) = 1/sqrt(2)|0&gt; - 1/sqrt(2)|1&gt;


So our transformation function is simply:

f(|00&gt;) = (1/sqrt(2)|0&gt; + 1/sqrt(2)|1&gt;) ⊗ |0&gt;
f(|01&gt;) = (1/sqrt(2)|0&gt; + 1/sqrt(2)|1&gt;) ⊗ |1&gt;
f(|10&gt;) = (1/sqrt(2)|0&gt; - 1/sqrt(2)|1&gt;) ⊗ |0&gt;
f(|11&gt;) = (1/sqrt(2)|0&gt; - 1/sqrt(2)|1&gt;) ⊗ |1&gt;


Expand this out to our normalized probability vector form...

f(|00&gt;) = [ 1/sqrt(2), 1/sqrt(2) ] ⊗ [ 1, 0 ]
f(|01&gt;) = [ 1/sqrt(2), 1/sqrt(2) ] ⊗ [ 0, 1 ]
f(|10&gt;) = [ 1/sqrt(2), -1/sqrt(2) ] ⊗ [ 1, 0 ]
f(|11&gt;) = [ 1/sqrt(2), -1/sqrt(2) ] ⊗ [ 0, 1 ]


Now let's actually solve this...

f(|00&gt;) = [ 1/sqrt(2), 0, 1/sqrt(2), 0 ]
f(|01&gt;) = [ 0, 1/sqrt(2), 0, 1/sqrt(2) ]
f(|10&gt;) = [ 1/sqrt(2), 0, -1/sqrt(2), 0 ]
f(|11&gt;) = [ 0, 1/sqrt(2), 0, -1/sqrt(2) ]


That's our transformation function. 

Our  logic gate matrix, "lgm", is of size 2^n by 2^n where n = number of qubits being simulated, so in this case it's 2^2 by 2^2 or 4x4.
Our algorithm tells us that for each column i, set the column equal to f(i). We've defined our probability transformation function, so we can easily fill out these columns.

lgm = 
    |00&gt;       |01&gt;        |10&gt;        |11&gt;
[ 1/sqrt(2),         0,  1/sqrt(2),          0 ] 
[         0, 1/sqrt(2),          0,  1/sqrt(2) ]
[ 1/sqrt(2),         0, -1/sqrt(2),          0 ]
[         0, 1/sqrt(2),          0, -1/sqrt(2) ]


Now the final step in our algorithm is simply matrix-multiplying our matrix representing the entire quantum system, sys, by this logic gate, lgm.

And that does what we want. It will apply the hadamard gate only to the most qubit and leave the least significant qubit alone.
If you don't believe me, you can try it yourself and see that it works.

The reason this is so powerful is because it applies to any case.

Let's try this algorithm on your problem.

Imagine we have a 3 qubit system and we want to apply a CNOT gate to qubit[0] and qubit[2]. If you look up the CNOT matrix on Wikipedia, that matrix only applies to a 2-qubit system. A naive solution would be to append identity matrices to it using the Kronecker product to make it work on systems with three qubits. But this fails here: qubit[0] and qubit[2] are not adjacent, so simply appending identity matrices won't work.

We could swap qubit[0] with qubit[1], apply the CNOT gate, then swap them back. But that's slow. Instead, we could simply generate a non-adjacent CNOT gate for our problem using the algorithm above.

We first need to come up with a transformation function for each case. 

f(|000&gt;) = |0&gt; ⊗ |0&gt; ⊗ |0&gt;
f(|001&gt;) = |0&gt; ⊗ |0&gt; ⊗ |1&gt;
f(|010&gt;) = |0&gt; ⊗ |1&gt; ⊗ |0&gt;
f(|011&gt;) = |0&gt; ⊗ |1&gt; ⊗ |1&gt;
f(|100&gt;) = |1&gt; ⊗ |0&gt; ⊗ |1&gt;
f(|101&gt;) = |1&gt; ⊗ |0&gt; ⊗ |0&gt;
f(|110&gt;) = |1&gt; ⊗ |1&gt; ⊗ |1&gt;
f(|111&gt;) = |1&gt; ⊗ |1&gt; ⊗ |0&gt;


If you understand the CNOT gate, you can understand why this is our function. Think of this like a truth table. Since our control qubit is the most significant qubit, qubit[2], only when that qubit is |1> will the least significant qubit, qubit[0], be negated. 

Expand this out to our normalized probability vector form...

f(|000&gt;) = [ 1, 0 ] ⊗ [ 1, 0 ] ⊗ [ 1, 0 ]
f(|001&gt;) = [ 1, 0 ] ⊗ [ 1, 0 ] ⊗ [ 0, 1 ]
f(|010&gt;) = [ 1, 0 ] ⊗ [ 0, 1 ] ⊗ [ 1, 0 ]
f(|011&gt;) = [ 1, 0 ] ⊗ [ 0, 1 ] ⊗ [ 0, 1 ]
f(|100&gt;) = [ 0, 1 ] ⊗ [ 1, 0 ] ⊗ [ 0, 1 ]
f(|101&gt;) = [ 0, 1 ] ⊗ [ 1, 0 ] ⊗ [ 1, 0 ]
f(|110&gt;) = [ 0, 1 ] ⊗ [ 0, 1 ] ⊗ [ 0, 1 ]
f(|111&gt;) = [ 0, 1 ] ⊗ [ 0, 1 ] ⊗ [ 1, 0 ]


Now let's actually solve this...

f(|000&gt;) = [ 1, 0, 0, 0, 0, 0, 0, 0 ]
f(|001&gt;) = [ 0, 1, 0, 0, 0, 0, 0, 0 ]
f(|010&gt;) = [ 0, 0, 1, 0, 0, 0, 0, 0 ]
f(|011&gt;) = [ 0, 0, 0, 1, 0, 0, 0, 0 ]
f(|100&gt;) = [ 0, 0, 0, 0, 0, 1, 0, 0 ]
f(|101&gt;) = [ 0, 0, 0, 0, 1, 0, 0, 0 ]
f(|110&gt;) = [ 0, 0, 0, 0, 0, 0, 0, 1 ]
f(|111&gt;) = [ 0, 0, 0, 0, 0, 0, 1, 0 ]


Let's make these the columns of our lgm matrix.

lgm =
[ 1, 0, 0, 0, 0, 0, 0, 0 ]
[ 0, 1, 0, 0, 0, 0, 0, 0 ]
[ 0, 0, 1, 0, 0, 0, 0, 0 ]
[ 0, 0, 0, 1, 0, 0, 0, 0 ]
[ 0, 0, 0, 0, 0, 1, 0, 0 ]
[ 0, 0, 0, 0, 1, 0, 0, 0 ]
[ 0, 0, 0, 0, 0, 0, 0, 1 ]
[ 0, 0, 0, 0, 0, 0, 1, 0 ]


Now if the matrix-multiply our entire system matrix, sys, by our logic gate matrix, lgm, our result will be the effect of applying the CNOT gate to qubit[2] and qubit[0] where qubit[2] is the control qubit.