boundaries = getresult('HEAT','boundaries');

I = boundaries.I_bias_right;
V = boundaries.V_bias_left;

plot(V,-I,'Voltage (V)','Current (A)');

P1 = pinch(boundaries.P_lead_left);
P2 = pinch(boundaries.P_lead_right);
P3 = pinch(boundaries.P_substrate);
P4 = pinch(boundaries.P_top_surface);

Ptotal = abs(P1+P2+P3+P4);

R = 50; # ohm
Panalytic = V^2/R;

plot(V,Ptotal,Panalytic,'Voltage (V)','Power (W)');
legend('simulation','analytic');