diff --git a/lib/Graphical_load_flow.m b/lib/Graphical_load_flow.m
new file mode 100644
index 00000000..4cf73d9a
--- /dev/null
+++ b/lib/Graphical_load_flow.m
@@ -0,0 +1,65 @@
+%% GRAPHICAL REPRESENTATION OF POWER FLOW (DRAFT CODE)
+%The following code, written as a function, would take the power flow
+%solution as an input and it would output the plot with the graphical
+%representation.
+function Graphical_load_flow(results)
+%results is a structure containing the PF solution as obtained from
+%Matpower.
+
+%Define the number of buses of the system
+num_nodi=length(results.bus(:,1));
+
+%Define the source and the sink buses for each line
+s=results.branch(:,1)
+t=results.branch(:,2)
+
+%Define the real power and reactive power flows over each line
+P=results.branch(:,14)
+Q=results.branch(:,15)
+
+%Define the buses to which the generators are connected
+Gen_Bus=results.gen(:,1)
+
+%Define the power generated by each generator
+Gen_Power=results.gen(:,2)
+
+%Define the buses to which the loads are connected
+Load_Bus=results.bus(find(results.bus(:,3)~=0),1)
+
+%Define the power absorbed by each load
+Load_Power=results.bus(:,3)
+
+%This code is required for the graph to display the correct directions 
+%of power flows.
+Neg=find(P<0)
+[t(Neg) s(Neg)]=deal(s(Neg), t(Neg))
+
+%Create the graph of the power flow
+G=digraph(s,t,abs(P))
+
+%Create the plot of the graph with various options (they could be set by
+%the user as an input to the function)
+Plot=plot(G,'LineWidth',abs(G.Edges.Weight)*0.02,'Layout','force','NodeLabel',[1:num_nodi],'NodeFontSize',8,'EdgeLabel',abs(G.Edges.Weight))
+
+%Set the color of the branches proportional to the power flow
+colormap jet;
+Plot.EdgeCData=abs(G.Edges.Weight);
+colorbar
+
+%Change the colors of the buses based on the characteristics of the bus
+highlight(Plot,Gen_Bus,'NodeColor','r')
+hold on
+highlight(Plot,Load_Bus,'NodeColor','g')
+%Gen+Load
+GL_Bus=intersect(Gen_Bus,Load_Bus);
+if ~isempty(GL_Bus)
+    highlight(Plot,GL_Bus,'NodeColor','m')
+end
+
+%Increase the size of the bus based on the power generated/absorbed
+for i=1:length(Gen_Bus)
+    highlight(Plot,Gen_Bus(i),'MarkerSize',Gen_Power(Gen_Bus(i))/10)
+end
+for i=1:length(Load_Bus)
+    highlight(Plot,Load_Bus(i),'MarkerSize',Load_Power(Load_Bus(i))/10)
+end
\ No newline at end of file