A numerical model of neurovascular coupling (NVC) is presented based on neuronal activity coupled to vasodilation/contraction models via the astrocytic mediated perivascular K+ and the smooth muscle cell (SMC) Ca2+ pathway termed a neurovascular unit (NVU). Luminal agonists acting on P2Y receptors on the endothelial cell (EC) surface provide a flux of inositol trisphosphate (IP3) into the endothelial cytosol. This concentration of IP3 is transported via gap junctions between EC and SMC providing a 10 source of sacroplasmic derived Ca2+ in the SMC. The model is able to relate a neuronal input signal to the corresponding vessel reaction (contraction or dilation). A tissue slice consisting of blocks, each of which contain an NVU is connected to a space filling H-tree, simulating a perfusing arterial tree (vasculature) The model couples the NVUs to the vascular tree via a stretch mediated Ca2+ channel on both the EC and SMC. The SMC is induced to oscillate by increasing an agonist flux in the EC 15 and hence increased IP3 induced Ca2+ from the SMC stores with the resulting calcium-induced calcium release (CICR) oscillation inhibiting NVC thereby relating blood flow to vessel contraction and dilation following neuronal activation. The coupling between the vasculature and the set of NVUs is relatively weak for the case with agonist induced where only the Ca2+ in cells inside the activated area becomes oscillatory however, the radii of vessels both inside and outside the activated area oscillate (albeit small 20 for those outside). In addition the oscillation profile differs between coupled and decoupled states with the time required to refill the cytosol with decreasing Ca2+ and increasing frequency with coupling. The solution algorithm is shown to have excellent weak and strong scaling. Results have been generated for tissue slices containing up to 4096 blocks.