% Konstruktionsuppgift 1, uppg 5
clear all

Px=0.01;
my=1.8e-5;
h=0.06;

% Analytisk lösn
y=(0:0.03/100:0.03);
Uana=16.66*y-Px/(2*my)*y.^2;

for n=1:6
N=10^n;
dh=h/(2*N);
e1=ones(N,1);
e2=-2*ones(N,1);
e2(1)=1;
e2(N)=-1;
A=spdiags([e1 e2 e1],[-1 0 1],N,N);
b=ones([N,1]);
b=(-Px*dh^2/my).*b;
b(1)=0;
b(N)=0;
eval(['U' num2str(N) '=A\b;']);
eval(['y' num2str(N) '=(0:0.03/(length(U' num2str(N) ')-1):0.03);']);
eval(['e(n)=(Uana(length(Uana))-U' num2str(N) '(length(U'...
    num2str(N) ')))/U' num2str(N) '(length(U' num2str(N) '));']);
grid(n)=N;
end

subplot(2,1,1)
plot(U10,y10,'b.')
hold on
plot(U100,y100,'r.')
plot(U1000,y1000,'g')
plot(Uana,y,'k')
legend('10 steg','100 steg','1000 steg','analytisk lösning')
xlabel('hastighet [m/s]');
ylabel('Y - koordinat');

subplot(2,1,2)
loglog(grid,abs(e),'-*')
xlabel('antal indelningar');
ylabel('e - relativa felet');