To make Figs. 7 & 10: run the relevant parts of the script plotFigs7&10.m (see code for details). But that's just plotting the numbers. The important question is: how do you GET the numbers? We will go through that now. It's a complicated process, so be prepared...

First of all, you'll need the specific workspaces I have used. The numbers in Fig. 7 from H with W_E have been obtained from the data in:
SimCity/L50W50ff0_02: noise realisations 1, 2, 3, 4, 5
The numbers in Fig. 7 from H with V have been obtained from the same data, plus:
SimCity/L50W50ff0_1: noise realisation 1
The same exact data was used for the top and bottom panel. I've uploaded both the original networks (after running ValleyNetwork.m) and the processed workspaces (after running ProcessNetwork.m on the 'original' workspaces; these are in folders with names ending in 'P').

You'll need to know which eigenstates were used in each noise realisation, and the domain pairs examined for decay in each eigenstate. I will write n->m to indicate that we measured the decay coefficient from domain n to m. You can run DomLabels.m to plot and label domain numbers for visualisation and ease of use.

SimCity/L50W50ff0_02, noise realisation 1. evec 1: 1->2
SimCity/L50W50ff0_02, noise realisation 2. evec 1: 2->1, 2->3, 2->4
SimCity/L50W50ff0_02, noise realisation 3. evec 1: 1->3, 1->2
SimCity/L50W50ff0_02, noise realisation 4. evec 1: 1->3, 1->2. evec 2: 3->4, 3->2
SimCity/L50W50ff0_02, noise realisation 5. evec 1: 2->1, 2->3

SimCity/L50W50ff0_1, noise realisation 1.
evec 1: 12->11, 12->14, 12->19, 12->9, 12->6, 12->16, 14->18, 18->15, 11->8, 11->5
evec 2: 18->15, 18->13, 18->14, 14->11, 11->8, 11->5, 12->19, 12->9, 12->6
evec 3: 11->8, 11->5, 5->3, 11->6, 11->13, 11->14, 12->9, 12->19
evec 4: 19->16, 12->14, 12->9, 12->6, 11->8, 11->5, 19->12, 12->11
evec 5: 4->7, 4->3, 4->8, 7->10, 3->5

For Fig. 10, I've used data uploaded to SimCity/L50W50ff0_06: noise realisation 1. The specific nearest neighbour pairs used in the first few eigenstates are:

evec 2: 8->4, 7->4, 4->2
evec 3: 2->4
evec 4: 3->2, 2->4, 4->8
evec 5: 4->2, 4->8, 8->7, 2->3
evec 6: 2->3, 8->7, 4->7
evec 7: 2->3, 2->5, 4->7, 8->7

In general, OBVIOUSLY you don't have to use these specific examples. You should be able to run the test described below on ANY correctly selected domain pairs, in any eigenstate, in any noise realisation, and the conclusion should be unchanged (see discussion in the paper for how to choose suitable examples). The information given above is necessary only so that you can reproduce the exact numbers I've plotted.

OK, so now to measure the decay rate predicted by LLT. Let's pull up one of the relevant processed workspaces listed above as an example (i.e. load it into Matlab, please). There is a matrix called 'rhoij'. This contains the *MEAN* Agmon distance between different domains in the network at different energies. This should be used for generating the top panel of Fig. 7, and Fig. 10. For the bottom panel of Fig. 7, you have to do this: Open the code ProcessNetwork.m. Look at the section on lines 1565-1584 (labelled 'step 7' in the comments). This is where this rhoij matrix is calculated. At this point we want the approximation of the *TRUE* Agmon distance. Line 1580 reads:
rhoij(n,k) = mean(options);
If you change it to:
rhoij(n,k) = min(options);
this segment will calculate the what we want. So copy this segment of code, make this change, paste it into the command window and execute. Great. Now open (from the workspace) a matrix called 'nn'. This is a list of nearest neighbour domains. Note the row numbers of the domain pairs of interest. Next open a vector called 'Evec'. This is the energy vector in units of E_0 at which we have evaluated the Agmon distance. The matrix rhoij has rows corresponding to domain pairs in 'nn', and columns corresponding to the entries in 'Evec'. All you have to do now is read the values out for the correct domain pairs at the correct energies. The right energies are the eigenvalues of the corresponding eigenstates, which we'll talk about below. If the eigenvalue is not close to an entry of 'Evec' (e.g. it's in between), I used linear interpolation using the nearest two Evec entries to estimate rho.

So far, so good. Now we need to extract the same numbers from the eigenstates. We use the script RhoFromEsts.m. Detailed instructions on how to use it are included in the comments therein. Once you run it, you'll get the eigenvalue of the eigenvector you're inspecting, and a vector called 'weight' which contains the mean (average) of |\psi| on each domain. In this vector, the domains are ordered as per labels on the pictures generated by this code (also by DomLabels.m); these are also the domain labels quoted above. Say you want to know the decay coefficient from domain 's' to 't'. Then calc log(weight(s)/weight(t)), which will be equivalent to \rho (see eqn (8) in the paper).