%% radiale Laplace-Gleichung
%

% Daten
ri = 0.1; %Radius innen und extern
re = 10;
ui = -1; % innerer Randwert
ue = 0; % aussen-Randwert

n = 100; % innere Diskr.Punkte
h = (re-ri)/(n+1);

j= 1:n; %Indexvektor
rj = ri + h*j;

% Matrixdiagonalen
d = 2*ones(n,1);
d_oben = -(1+ h./(2*rj(1:end-1)));
d_unten = -(1- h./(2*rj(2:end)));

A = diag(d) + diag(d_oben,1) + diag(d_unten,-1);

% rechte Seite enthaelt Randterme
b = zeros(n,1);
b(1) = (1-h/(2*rj(1)))*ui;
b(n) = (1+h/(2*rj(n)))*ue;

% Loesung
u = A\b;

% Exakte Loesung
uex = (ui*log(rj/re) + ue*log(ri./rj))/log(ri/re);

% Plot

plot([ri,rj,re],[ui, u', ue], [ri,rj,re],[ui, uex, ue])

figure(2)
plot([ri,rj,re],[0, u'-uex, 0] )


