
T0 = 1;
kappa = 1/12;
tmax = 0.2;

ns = 10 * 2.^[0:5];

grids = length(ns);
e_a = zeros(1, grids);
e_b = zeros(1, grids);
e_c = zeros(1, grids);

for k = 1:length(ns)
  n = ns(k);
  disp(sprintf('n = %d', n))
  xs = [0:1/n:1];
  Tinit = [T0, zeros(1, n - 1), T0];
  lsg = waermeleitung_loesung(xs, tmax, T0, kappa);

  s = 1/24;
  l_a = wleit2(Tinit, s, round(tmax * kappa / s / (1/n)^2));
  e_a(k) = l1diff(l_a, lsg);

  s = 1/6;
  l_b = wleit2(Tinit, s, round(tmax * kappa / s / (1/n)^2));
  e_b(k) = l1diff(l_b, lsg);

  s = 5/12;
  l_c = wleit2(Tinit, s, round(tmax * kappa / s / (1/n)^2));
  e_c(k) = l1diff(l_c, lsg);
end

quot_a = e_a(1:grids-1) ./ e_a(2:grids);
order_a = log(quot_a) / log(2);

quot_b = e_b(1:grids-1) ./ e_b(2:grids);
order_b = log(quot_b) / log(2);

quot_c = e_c(1:grids-1) ./ e_c(2:grids);
order_c = log(quot_c) / log(2);

disp('** s = 1/24: Fehler / Quotient / Ordnung')
disp(sprintf(' %g', e_a))
disp(quot_a)
disp(order_a)

disp('** s = 1/6:  Fehler / Quotient / Ordnung')
disp(sprintf(' %g', e_b))
disp(quot_b)
disp(order_b)

disp('** s = 5/12: Fehler / Quotient / Ordnung')
disp(sprintf(' %g', e_c))
disp(quot_c)
disp(order_c)

figure(1)
loglog(ns, e_a, '-*', ns, e_b, '-o', ns, e_c, '-s')
legend('s = 1/24', 's = 1/6', 's = 5/12')
xlabel('Gitterpunkte')
ylabel('L_1 Fehler')

%figure(2)
%plot(ns(2:grids), order_a, ns(2:grids), order_b, ns(2:grids), order_c)
%legend('s = 1/24', 's = 1/6', 's = 5/12')
%xlabel('Gitterpunkte')
%ylabel('Ordnung')
