I fail it for the plot.

% Metody teoretyczne

clear
N=1; T=2; % falsz, prawda

%zadawanie wartosci prawdopodobienstw

%wirus
pw(T)=0.001;
pw(N)=1-pw(T);

% błąd użytkownika
pu(T)=0.5;
pu(N)=1-pu(T);

% wysoka temp
pt(T)=0.05;
pt(N)=1-pt(T);

% błąd sys
ps(T,N,N)=0.005;
ps(T,N,T)=0.007;
ps(T,T,N)=0.55;
ps(T,T,T)=0.6;
ps(N,N,N)=1-ps(T,N,N);
ps(N,N,T)=1-ps(T,N,T);
ps(N,T,N)=1-ps(T,T,N);
ps(N,T,T)=1-ps(T,T,T);

% awaria sprzętu
pa(T,T)=0.1;
pa(T,N)=0.01;
pa(N,T)=1-pa(T,T);
pa(N,N)=1-pa(T,N);

%zawieszenie
pz(T,T,T)=0.95;
pz(T,N,T)=0.90;
pz(T,T,N)=0.40;
pz(T,N,N)=0.001;
pz(N,T,T)=1-pz(T,T,T);
pz(N,N,T)=1-pz(T,N,T);
pz(N,T,N)=1-pz(T,T,N);
pz(N,N,N)=1-pz(T,N,N);

% tablica prawdopodobienstw lacznych

for S=1:2
  for A=1:2
    for Z=1:2
      for U=1:2
        for W=1:2
            for T=1:2
               p(S,A,Z,U,W,T)=pw(W)*pu(U)*pt(T)*ps(S,W,U)*pa(A,T)*pz(Z,S,A);
            end
        end
      end
    end
  end
end

% Sprawdzenie czy suma wszystkich prawdopodobienstw jest rowna 1
sum(sum(sum(sum(sum(sum(p))))))


% Na podstawie tablicy prawdopodopienstw lacznych mozna obliczyc dowolne prawdopodobienstwo, np.P(A|S)=P(A,S)/P(S)

PA=sum(sum(sum(sum(sum(p(:,2,:,:,:,:))))))
PUZ=sum(sum(sum(sum(p(:,:,2,2,:,:)))))/sum(sum(sum(sum(sum(p(:,:,2,:,:,:))))))

% -------------------------------------------------------
% Metody Monte Carlo

K=10000; % jednorazowa porcja taktow
close
hold on                             % polecenie zostawia wykres, nastepne beda rysowane na tym samym
plot([1 100],[PUZ PUZ])          % wykres P(A|S) wyznaczonego teoretycznie
axis([0 100 0 1])                % zadanie skali rysunku

PLOCIK=[]
SR=[]

for j=1:100
    w=rand(1,K)<pw(T);
    u=rand(1,K)<pu(T);
    sNN=rand(1,K)<ps(T,N,N);
    sNT=rand(1,K)<ps(T,N,T);
    sTN=rand(1,K)<ps(T,T,N);
    sTT=rand(1,K)<ps(T,T,T);
    s=(w&u&sTT)|(w&~u&sTN)|(~w&u&sNT)|(~w&~u&sNN);

    t=rand(1,K)<pt(T);

    aN=rand(1,K)<pa(T,N);
    aT=rand(1,K)<pa(T,T);
    a=(t&aT)|(~t&aN);

    zNN=rand(1,K)<pz(T,N,N);
    zNT=rand(1,K)<pz(T,N,T);
    zTN=rand(1,K)<pz(T,T,N);
    zTT=rand(1,K)<pz(T,T,T);
    z=(s&a&zTT)|(s&~a&zTN)|(~s&a&zNT)|(~s&~a&zNN);

    temp=sum(u&z)/sum(z)
    PLOCIK(j)=temp
    temp2=0;
    for i=1:length(PLOCIK)
        temp2=temp2+PLOCIK(i);
    end
    temp2=temp2/length(PLOCIK);
    SR(j)=temp2;
    % obliczenie czestosci w j-tym kroku na podstawie odpowiednich tablic jedynek, np. P(A|S)= sum(a&s)/sum(s)

    % usrednienie czestosci dla wszystkich dotychczasowych krokow

    % narysowanie wykresu: plot([1:j], tablica_srednich_czestosci)


end

    plot([1:j], SR(1:j))