clc;
format long;
a=10^(-3);
Us=4;
g=10^(-4);
R=10^4;
E0=75;

%graficznie

u=linspace(0,80,1000);
i=a*(u.^2).*exp(-u/Us)+g*u;
f=E0-R*i;
plot(i,u,i,f)
hold on;
grid on;
xlabel('I [A]')
ylabel('U [V]')
title('Punkty pracy rezystora nieliniowego - metoda graficzna');

%numerycznie - metoda Nortona
clear u;
clear i;
u=linspace(0,80,1000);
%rownanie wynikajace z rozwiazania zadania metoda Nortona
y=-E0./R+u./R+a.*u.^2.*exp(-u./Us)+g.*u;
figure(2)
y1=u-u;
axis([-0.008 0.010 0 80])
plot(u,y,u,y1)
%widac, ze pierwiastki sa w okolicach 5,15 i 35
%wokol tych punktow szukamym miejsc zerowych funkcja fzero
ylabel('10U^2exp(-U/4)+2U-75')
xlabel('U')
title('Wykres wynikajacy z rozwiazania ukladu metoda Nortona')
hold on;
grid on;
fun=inline('u.*u.*10.*exp(-u./4)+2.*u-75');
%napiecia
U1=fzero(fun, 10)
U2=fzero(fun, 20)
U3=fzero(fun, 40)
%i odpowiadajace im prady w punkcie pracy
I1=1000*(a*U1.^2.*exp(-U1/Us)+g*U1)
I2=1000*(a*U2.^2.*exp(-U2/Us)+g*U2)
I3=1000*(a*U3.^2.*exp(-U3/Us)+g*U3)