% Richtungsfeld einer dgl

% Länge des Steigungsvektors
l=0.5;
%y-bereich
Y=-10:0.5:10;
% x-bereich
X=-10:0.5:10;
%dgl y'=f(y,x)
dgl = @(y, x) sqrt(1+y^2);



X1=2:.1:8;
X2=-3:0.1:3;
X3=-8:0.1:-2;
counter=0;

X_RESULT=[];
Y_RESULT=[];
tic;
x1 = zeros(1, numel(X)*numel(Y));
x2 = zeros(1, numel(X)*numel(Y));
y1 = zeros(1, numel(X)*numel(Y));
y2 = zeros(1, numel(X)*numel(Y));
for x0 = X
    for y0 = Y
        m=dgl(y0, x0);
        deltaX=l/sqrt(1+m^2);
        deltaY=m*deltaX;
        counter=counter+1;
        x1(:, counter) =x0-deltaX/2;
        x2(:, counter) =x0+deltaX/2;
        y1(:, counter) =y0-deltaY/2;
        y2(:, counter) =y0+deltaY/2;
    end
end
plot([x1; x2], [y1; y2], 'b');
hold on
toc;
plot( X2, sinh(X2), X1, sinh(X1-5), X3, sinh(X3+5));
hold off
xlim([-10,10]);
grid on;

ylim([min(Y),max(Y)]);


counter