clc;
fprintf('1. Rysowanie wykresu pochodnej funkcji f(x) na przedziale [x1,x2]\n');
fprintf('2. Wyznaczanie pochodnej funkcji f(x) w punkcie x0\n');
fprintf('3. Rysowanie wykresow pochodnych czastkowych funkcji f(x,y) na przedziale [x1,x2]x[y1,y2]\n');
fprintf('4. Wyznaczanie gradientu funkcji f(x,y) w punkcie (x0,y0)\n');
fprintf('5. Wyznaczanie gradientu funkcji f(x,y,z) w punkcie (x0,y0,z0)\n\n');
k=input('Wybierz (1-5): ');
h=0.001; % przyrost w def. pochodnej dla case 2,4,5
switch k
  case 1
    x1=input('x1=');
    x2=input('x2=');
    x=linspace(x1,x2);
    f=input('f(x)=');
    if (length(f)==1) % jesli f jest funkcja stala
      f=f*ones(1,100);
    end
    Df=diff(f)./diff(x);
    plot(x(1:99),[f(1:99);Df]);
    xlabel('x');
    ylabel('y');
    legend('f(x)','Df(x)');
  case 2
    x0=input('x0=');
    x=[x0-h/2,x0+h/2];
    f=input('f(x)=');
    if (length(f)==1) % jesli f jest funkcja stala
      f=f*ones(1,2);
    end
    Df=diff(f)/h;
    fprintf('Df(%f)=%f\n',x0,Df);
  case 3
    x1=input('x1=');
    x2=input('x2=');
    y1=input('y1=');
    y2=input('y2=');
    x=linspace(x1,x2);
    y=linspace(y1,y2);
    [x,y]=meshgrid(x,y);
    f=input('f(x,y)=');
    if (length(f)==1) % jesli f jest funkcja stala
      f=f*ones(100,100);
    end
    Dfx=transpose(diff(transpose(f))./diff(transpose(x)));
    Dfy=diff(f)./diff(y);
    subplot(1,3,1); % wykres f(x,y)
    plot3(x,y,f);
    title('f(x,y)');
    xlabel('x');
    ylabel('y');
    zlabel('z');
    subplot(1,3,2); % wykres df/dx(x,y)
    plot3(x(:,1:99),y(:,1:99),Dfx);
    title('df/dx(x,y)');
    xlabel('x');
    ylabel('y');
    zlabel('z');
    subplot(1,3,3); % wykres df/dy(x,y)
    plot3(x(1:99,:),y(1:99,:),Dfy);
    title('df/dy(x,y)');
    xlabel('x');
    ylabel('y');
    zlabel('z');
  case 4
    x0=input('x0=');
    y0=input('y0=');
    x=[x0-h/2,x0,x0+h/2];
    y=[y0-h/2,y0,y0+h/2];
    [x,y]=meshgrid(x,y);
    f=input('f(x,y)=');
    if (length(f)==1) % jesli f jest funkcja stala
      f=f*ones(3,3);
    end
    Dfx=(f(2,3)-f(2,1))/h; % tzn. (f(x0+h/2,y0)-f(x0-h/2,y0))/h
    Dfy=(f(3,2)-f(1,2))/h; % tzn. (f(x0,y0+h/2)-f(x0,y0-h/2))/h
    Df=[Dfx,Dfy];
    fprintf('Df(%f,%f)=(%f,%f)\n',x0,y0,Dfx,Dfy);
  case 5
    x0=input('x0=');
    y0=input('y0=');
    z0=input('z0=');
    x=[x0-h/2,x0,x0+h/2];
    y=[y0-h/2,y0,y0+h/2];
    z=[z0-h/2,z0,z0+h/2];
    [x,y,z]=meshgrid(x,y,z);
    f=input('f(x,y,z)=');
    if (length(f)==1) % jesli f jest funkcja stala
      f=f*ones(3,3,3);
    end
    Dfx=(f(2,3,2)-f(2,1,2))/h; % tzn. (f(x0+h/2,y0,z0)-f(x0-h/2,y0,z0))/h
    Dfy=(f(3,2,2)-f(1,2,2))/h; % tzn. (f(x0,y0+h/2,z0)-f(x0,y0-h/2,z0))/h
    Dfz=(f(2,2,3)-f(2,2,1))/h; % tzn. (f(x0,y0,z0+h/2)-f(x0,y0,z0-h/2))/h
    Df=[Dfx,Dfy,Dfz];
    fprintf('Df(%f,%f,%f)=(%f,%f,%f)\n',x0,y0,z0,Dfx,Dfy,Dfz);
  otherwise
    fprintf('Blednie wybrane polecenie\n');
end % koniec switcha

% eof