clear;
close all;
skellam = Skellam(1024,0,32);
k=1;
Q=20;
xq=linspace(-Q,Q,1024);
for h=4:0.01:10
    j=1;
    for x=-20:0.6:20;
        suma=0;
        for i=1:length(skellam)
            if ((((skellam(i)/h)-floor(x/h))>=0) && (((skellam(i)/h)-floor(x/h))<1))
               suma=suma+1/h;
            end
        end
        wynik1(k,j)=suma/length(skellam);
        j=j+1;
    end
    j=1;
    for x=-20:0.6:20;
        suma=0;
        for i=1:length(skellam)
            if ((((skellam(i)/h)-(x/h))>=0) && (((skellam(i)/h)-(x/h))<1))
               suma=suma+1/h;
            end
        end
        wynik2(k,j)=suma/length(skellam);
        j=j+1;
    end
    suma1=0;
    for l=1:length(xq)
        suma1 = suma1+(wynik1(k,floor(xq(i)))-f(xq(i)))^2;
    end
    val1(k)=suma1/(2*Q);
    suma2=0;
    for l=1:length(xq)
        suma2 = suma2+(wynik2(k,floor(xq(i)))-f(xq(i)))^2;
    end
    val2(k)=suma2/(2*Q);
    k=k+1;
end
[a1 b1] = min(val1);
[a1 b2] = min(val2);
h=4:0.01:10;
h(b2)
h(b1)