nl = 10;
n = 10;
W = rand(n,n,nl);
b = rand(n, nl);
xv = rand(n,1);
yv = rand(n,1);

dW = zeros(size(W));
db = zeros(size(b));

[z,y,z_0,val0] = f_forward(W, b, xv, yv);
temp = 2*(z_0-yv).*sigmader(y(:,1));
for k=1:(nl-1)
    dW(:,:,k) = temp*z(:,k)'; 
    db(:,k) = temp;
    temp = (temp'*W(:,:,k).*sigmader(y(:,k+1))')';
end
dW(:,:,nl) = temp*z(:,nl)'; db(:,nl) = temp;

[z,y,z_0,val] = f_forward(W-dW, b-db, xv, yv);
val0
val

function y=sigma(x)
    y = 1./(1+exp(-x));
end

function y = sigmader(x)
    y = sigma(x).*(1-sigma(x));
end