My expectation is that the second part is just requiring you to go beyond only using change in bus voltage magnitude at each step (to determine when the case is solved) and adding a check for the \$P\$ and \$Q\$ tolerance. I would think you would want to calculate the total \$P\$ and \$Q\$ at your bus 2 and sum it with the load injection at that bus. If you have a perfect solution the result will be zero. If not, iterate again.

Using this figure as a reference,

The power injection at your bus i will be,

$$P=\sum_{j=1}^n|V_i||Y_{ij}||V_j|cos(\delta _j+\theta_{ij}-\delta_i) $$

From this you subtract the scheduled load for that bus, \$P_i\$, and the closer to zero the better the solution. This is what I think your \$\Delta P\$ will be.

Do similar for reactive with,

$$Q=\sum_{j=1}^n|V_i||Y_{ij}||V_j|sin(\delta _j+\theta_{ij}-\delta_i)$$

I highly recommend you use another tool to check yourself. MATPOWER is free and runs on top of Matlab (or Octave if you prefer a free equivalent to Matlab). I created a test case file like yours (I named it test.m) and set the bus voltage base to 1.0kV to leave it in per-unit.

function mpc = test
%% MATPOWER Case Format : Version 2
mpc.version = '2';
%%-----  Power Flow Data  -----%%
%% system MVA base
% This is the base that the branch impedances are on.  
mpc.baseMVA = 100;
%% Vm under mpc.gen sets the slack bus (1) voltage magnitude
%% Va under mpc.bus sets the slack bus (1) voltage angle 
%   bus_i   type      Pd       Qd      Gs       Bs  area      Vm      Va      baseKV    zone     Vmax      Vmin
mpc.bus = [
  1    3     0      0     0      0     1     1     0      1.0     1      1       1;
  2    1  0.90   0.25     0      0     0     1     1      1.0     1      1       0;
];
%% generator data
%   bus    Pg       Qg    Qmax    Qmin        Vg        mBase   status  Pmax    Pmin    Pc1 Pc2 Qc1min  Qc1max  Qc2min  Qc2max  ramp_agc    ramp_10 ramp_30 ramp_q  apf   ?
mpc.gen = [ 
     1      0       0        0       0      1.05     0       1    0     0    0   0       0       0       0       0         0        0       0      0    0   0;
];
%% branch data
%   fbus    tbus             r           x      b      rateA    rateB   rateC    ratio   angle   status   angmin      angmax
mpc.branch = [
 1     2         0.0        0.12   0        0      0     0      1        0        1        0         0;
];

I ran it with the runpf ('test') command in Octave and the results from Newton Rahpson load flow was bus 2 voltage of 1.05∠-0.056° pu. I didn't try real hard to check the data, I just wanted to point you to a way to check your own power flow results.

Also, be careful with your code notes, could trip you up later. For example you have:

% conjugate ((P(2) + jQ(2)/VK(2)), so:

But I think you mean, "% conjugate [(P(2) + jQ(2)]/VK(2), so:"