Is my backpropagation code correct?
Thanks for everyone's help. I have now solved the problem. I hadn't quite understood the back propagation algorithm. I invite people to take a look at this link: backpropagation @ AGH UST which has solved my problem.

Code:  

clear
pkg load image
graphics_toolkit("gnuplot")
tic
sigmoid = @(z) 1./(1 + exp(-z));
sig_der = @(y) sigmoid(y).*(1-sigmoid(y));


function [cost, mid_layer, last_layer] = forward(w1,w2,data,sigmoid,i)
  mid_layer(:,1)=sum(w1.*data(1:2,i));
  mid_layer(:,2)=sigmoid(mid_layer(:,1));
  last_layer(:,1)=sum(mid_layer(:,2).*w2);
  last_layer(:,2)=sigmoid(last_layer(:,1));
  cost = data(3,i)-last_layer(:,2);
end

function [w1, w2] = backprop(w1,w2,mid_layer,last_layer,data,cost,sig_der,sigmoid,i)
  delta(1) = cost;
  delta(2:4) = cost .* w2;
  w2 = w2 + 0.99 .* delta(1) .* mid_layer(:,2) .* sig_der(last_layer(:,1));
  w1 = w1 + 0.99 .* delta(2:4) .* sig_der(mid_layer(:,1))' .* data(1:2,i);
  %b2 = b2 + 0.01 .* 
  %b1 = b1 + 0.01 .*
end

tic


data(:,1)=[0; 0; 0];
data(:,2)=[1; 0; 1];
data(:,3)=[0; 1; 1];
data(:,4)=[1; 1; 0];

w1=rand(2,3)./2.*(rand(2,3).*-2+1);
w2=rand(3,1)./2.*(rand(3,1).*-2+1);
%b1=ones(1,3);%rand(1,3)./2;
%b2=ones(1,1);%rand(1,1)./2;

for j=1:10000
  for i=[randperm(4)]
    [cost, mid_layer, last_layer] = forward(w1,w2,data,sigmoid,i);
    [w1, w2] = backprop(w1,w2,mid_layer,last_layer,data,cost,sig_der,sigmoid,i);
    cost_mem(j,i)=cost;
  end
end

toc

%{
dlmwrite("/tmp/weights_1",w1);
dlmwrite("/tmp/weights_2",w2);
%}