a=0; %wartość początkowa przedziału dla x
b=1; %wartość końcowa przedziału dla x
h1=0.1; %wartość kroku
h2=0.05;
h3=0.025;
N1=(b-a)/h1; %ilość iteracji
N2=(b-a)/h2;
N3=(b-a)/h3;
x1(1)=a; %pierwsza wartość z przedziału
x2(1)=a;
x3(1)=a;
y1(1)=1; %pierwsza wartość funkcji y (Runge-Kutty)
y2(1)=1;
y3(1)=1;
z1(1)=1; %pierwsza wartość funkcji z (Euler)
z2(1)=1;
z3(1)=1;
d1(1)=1; %pierwsza wartość funkcji d (dokładne)
d2(1)=1;
d3(1)=1;
for i=1:N1
    x1(i+1)=x1(1)+i*h1; %wartości wektora x
    z1(i+1)=z1(i)+h1*fun(x1(i),z1(i));  %rozwiązanie dla metody Eulera
    y1(i+1)=y1(i)+h1*fun(x1(i)+h1/2,y1(i)+h1/2*fun(x1(i),y1(i))); %rozwiązanie dla metody Rungego-Kutty
    d1(i+1)=2*exp(x1(i+1))-x1(i+1)-1; %rozwiązanie dokładne
end
for i=1:N2
    x2(i+1)=x2(1)+i*h2; %wartości wektora x
    z2(i+1)=z2(i)+h2*fun(x2(i),z2(i));  %rozwiązanie dla metody Eulera
    y2(i+1)=y2(i)+h2*fun(x2(i)+h2/2,y2(i)+h2/2*fun(x2(i),y2(i))); %rozwiązanie dla metody Rungego-Kutty
    d2(i+1)=2*exp(x2(i+1))-x2(i+1)-1; %rozwiązanie dokładne
end
for i=1:N3
    x3(i+1)=x3(1)+i*h3; %wartości wektora x
    z3(i+1)=z3(i)+h3*fun(x3(i),z3(i));  %rozwiązanie dla metody Eulera
    y3(i+1)=y3(i)+h3*fun(x3(i)+h3/2,y3(i)+h3/2*fun(x3(i),y3(i))); %rozwiązanie dla metody Rungego-Kutty
    d3(i+1)=2*exp(x3(i+1))-x3(i+1)-1; %rozwiązanie dokładne
end